Recreational and Nature-Based Tourism as a Cultural Ecosystem Service. Assessment and Mapping in a Rural-Urban Gradient of Central Spain

: Land management focused from the social-ecological perspective of ecosystem services should consider cultural services in decision-making processes. Nature-based tourism offers a great potential for landscape conservation, local development and the well-being of human populations. However, the subjectivity of recreational ecosystem services has meant a clear impediment to assessing and mapping them. In this study, an integrated numerical spatial method is developed, which quantiﬁes the supply and demand of recreational ecosystem services and allows mapping their spatial correspondence along a rural-urban gradient. The procedure also allows quantifying the inﬂuence of the landscape structure and the presence of protected areas on the degree of coupling between supply of recreational ecosystem services and demand for outdoor recreation and nature-based tourism and reveals that protected areas are hotspots of recreational ecosystem services. The results obtained highlight the usefulness of the methodological procedure developed as a tool for sustainable land planning and management from an integrative social-ecological approach.


Introduction
Cultural rural landscapes are characterised by their multi-functionality, stability and adaptation to local conditions [1]. Traditional agricultural activities have led to different types of land uses that have shaped complex landscape patterns with a high associated biocultural diversity and providing multiple ecosystem services (ES) [2]. The concept of ecosystem services links the functioning of ecosystems with social-ecological systems and human well-being [3,4]. Their supply is vulnerable to human use and therefore highly dependent on the functions of rural and urban ecosystems [5]. In this social-ecological context, the approach of rural-urban gradients is widely used to describe the relationships between land uses, urban sprawl, and ecological and socioeconomic dynamics [6][7][8] along with the transition from human-dominated landscapes to lightly developed rural areas [9,10]. These processes are the main driving forces of landscape structure change, impacting the supply and demand for ES. A comparison of ES supply in urban and rural regions can provide important arguments for establishing sustainable spatial planning strategies [11]. Therefore, a key challenge of land planning and management is the development of strategies based on the supply of multiple ES to control and regulate land use and urban development [12]. Likewise, the establishment of protected area networks (PAs), specially designed to preserve biodiversity and ecological flows [13], also have the potential to maintain or improve the supply of ES [14,15]. metropolitan area (Figure 1). This gradient is mainly characterized by pasture systems and forest-dominated landscapes in the northern uplands and agricultural systems in the southern lowlands. For several decades, the Madrid region has experienced notable socio-economic and landscape changes mainly related to intense processes of urban growth, urban sprawl [36] and rural abandonment, linked to the decrease of the population dedicated to agriculture and livestock and with the increase of employment in the secondary and tertiary sectors [34,37]. The widespread abandonment of agricultural and silvopastoral land uses and practices has favored a remarkable process of shrub encroachment and afforestation in rural areas [1,2,38]. These two main trends of landscape change, urban expansion and rurality loss, have caused an accentuated rural-urban dichotomy and significant social-ecological modifications in the study region [35]. However, in this transformed and dynamic territorial matrix, large areas with high natural and cultural values still persist, which have been recognized both regionally, nationally and internationally. Thus, a complex PA network that includes different protection categories has been established in the rural-urban gradient studied, covering one-third of the area ( Figure 1): (i) two Regional Parks ("Upper Manzanares River Basin Regional Park" and "Southeast Regional Park"), a protection status recognized by the Madrid Regional Government similar to the Category V of the International Union for Conservation of Nature (IUCN), aimed at protecting and conserving both natural values and those created by the interaction with human beings through traditional management practices; (ii) a National Park ("Sierra de Guadarrama National Park"), declared of general interest by the Government of Spain due to the high value of its well preserved natural systems; (iii) a Biosphere Reserve ("Upper Basins of the Manzanares, Lozoya and Guadarrama rivers"), an international category of protected area designated by UNESCO (Man and the Biosphere Program, MaB), which represents a balanced relationship between people and nature; (iv) six Nature 2000 sites (three Special Protection Areas, SPAs, designated under the Birds Directive of the European Union and three Sites of Community Importance, SCIs, determined by the Habitats Directive of the European Union). Both the protected landscapes of the study area and their surrounding lands have a great capacity to satisfy the recreational demand of NBT visitors [39].

Data Collection and Analyses
We use different types of data, collected both from public databases and available cartography and from questionnaires made to visitors to the area. As units of spatial analysis, we considered the 36 municipalities of the study area since they are the smallest governance unit in the Madrid Region and also the administrative scale of greater detail in which data are available from the socioeconomic and agricultural census (Refs. [1,2,34,35,40] among others).
The methodological approach was developed in accordance with the proposed objectives. Figure 2 summarizes the main steps followed, which are detailed below. The procedure followed for the assessment of the potential capacity of the land-scape to supply RES starts from the premise that the supply of ES is a function of the existing land use-land cover data (LULCs) in the territory [41]. According to several authors, we define "land cover" as the biophysical characteristics of a territory and "land use" as the utility that people obtain from it. Consequently, understanding the ways in which people use land is necessary to identify the potential supply of ES. The use of LULCs as proxies (or indirect indicators) of the ES supply is very frequent, due to the lack of empirical in-formation on ES flows [42]. Thus, proxy-based maps are more common and available that maps based on primary data [14,43]. In this work, valid and useful proxies were selected to evaluate the capacity of the landscape to supply RES. Proxies were land uses easily recognizable by landscape observers (such as agricultural systems, grasslands, forests or shrublands, among others). These landscape features can be considered as recreational resources and potential tourism attractors [23,39] and, therefore, as ES. Examples of landscape RES are forest systems, pastures, croplands and other natural and semi-natural habitats, which provide numerous social and cultural services and represent a privileged place for outdoor recreation and leisure [44,45]. Thus, we used 16 types of LULCs for the whole rural-urban gradient, derived from the reclassification of the CORINE Land Cover data set (2012) into more significant and representative land use categories, according to the dynamics of land uses characteristic of the studied region [46] and their capacity to provide opportunities for recreational activities and NBT [21] (Table 1). In each municipality of the gradient, we calculated the spatial coverage (%) of selected LULCs, considered as recreational services. We structured this information in a data matrix in which the municipalities were row vectors that contained quantitative data about the proxies of RES (column vectors). The dimensions of the final data matrix were 36 municipalities x 16 RES proxies ( Figure 3). Table 1. LULCs used as proxies for recreational and nature-based tourism ecosystem services (RES) supplied by the landscape. LULCs were quantified as a percentage of occupied area in each municipality of the rural-urban gradient. These proxies were also questions in a survey on landscape preferences of visitors to the gradient studied. In this survey, visitors qualitatively indicated which elements of the landscape they found attractive.
Savin juniper and Juniper forests 4.
Kermes oak and calcicolous shrublands 9. (d) each municipality of the rural-urban gradient is a row vector containing quantitative data (aij) of RES proxies (column data) in a matrix of recreational landscape services.

Step 2. Analysis of Outdoor Recreation Demand
LULCs, previously selected as RES proxies, were included as questions in a survey on landscape preferences of visitors to the study area. Thus, the assessment of landscape features by visitors allows to identify their RES demand [27]. The questionnaires consisted of a limited number of questions that could be easily and quickly answered and were conducted on weekends, holidays and vacation periods. The interviewees did not belong to any specific type of visitors; any visitor to the study area was considered in the survey. The respondents had to indicate which aspects of the landscape of those collected in the questionnaire (corresponding to the 16 RES selected; see Table 1) were most attractive to them, by marking the items that most closely matched their preferences. The structure and orientation of the questionnaires had been validated in previous works, with satisfactory results [22,39,[47][48][49].
The sampling covered the heterogeneity of the landscape of the studied gradient. It was carried out along routes and at rest points in selected areas considered to be of particular interest for tourism due to their natural, recreational and historical-cultural heritage values (towns and villages, tourist information points, visitor or interpretive centres of protected areas and picnic and camping zones). To avoid possible redundancies in the responses or overrepresentations in the responses, only one person from each group of visitors was randomly chosen to fill out the questionnaire [23,39]. In the sampling, a total of 400 contacts were made, of which 367 were considered to be valid for subsequent analyses.
Quantification of demand was based on the calculation of the frequency of positive responses to the questions of the survey. Thus, we obtained a column vector containing the score of the different items of the survey, i.e., quantitative information on the demand of visitors for the RES of the landscape.
Local people were excluded from the survey because they are the "insiders", the real managers of their landscape and its tangible and intangible attributes [50]. Therefore, they were not considered in the survey so that the sample was not biased by their possible interests, problems or overvaluations [23,51,52].

Step 3. Spatial Coupling between Supply and Demand for Landscape Recreational Services. Quantification and Mapping
We used the RES supply matrix (data aij; see Figures 3d and 4a) and the demand vector to calculate the correspondence between supply and demand for recreational services (data bij) along the rural-urban gradient ( Figure 4b). For this, we performed an algebraic product of matrix by vector (aij × bij). The resulting product vector informs about the degree of coupling supply-demand for RES and NBT at municipal level and its geospatial projection allowed us to obtain maps of this relationship (Figure 4c,d, respectively). The elements of this product vector (data cij) were divided into three categories (high, medium and low) by means of the natural break classification method ("Jenks optimization method" [53]), specifically designed to determine the best arrangement of values into different classes. Jenks's method is a within-groups variance minimization approach that determines class breaks to maximize homogeneity within classes and allow us to generate accurate maps, in terms of the representation of the spatial attributes of the data [54,55]. We thus obtained three types of municipalities along the rural-urban gradient, according to their greater or lesser correspondence between the supply and the demand for RES. We used ArcGis 10.6 (Esri, Redlands, CA, USA) as mapping tool.

Step 4. Analysis of the Relationship between Landscape Structure and the Supply and Demand for RES
We identified landscape structure by means of both landscape composition (using the selected LULCs, considered as valid proxies that captures the quantity of RES remaining in the landscapes; Ref. [41]) and landscape configuration (using landscape metrics, considered useful in providing an objective description of different aspects of landscape structure and patterns; Ref. [56]) Landscape metrics were selected according to criteria based on their comparability, non-redundancy, ease of interpretation, ability to describe spatial structure, and also on the satisfactory results provided by their use in previous studies [34,35,57]. The spatial metrics used were (Table 2): (1) Shannon's diversity index (SHDI, to quantify landscape diversity; it is a good indicator of landscape heterogeneity); (2) Patch richness index (PR, to measure the number of patches present); (3) Splitting index (SPLIT, to quantify landscape fragmentation [6]). Euclidean distance from the nearest neighbor (ENN, to describe the spatial isolation between patches and, consequently, the landscape connectivity). These metrics were calculated using the LULCs previously reclassified from Corine Land Cover map (2012). By means of a round moving window with a 100 m radius, we generated a raster map of each of these selected land metrics for the whole gradient studied. Subsequently, we extracted a mean value for each metric in each of the 36 municipalities included in the area. Metrics were calculated using Fragstats 4.2 [56] and were treated using ArcGis 10.6. Table 2. Landscape metrics used to calculate landscape spatial patterns. The calculation procedure and a description of each metric are indicated.

Landscape Metrics Formula Ranges Description
(1) Shannon's Diversity Index The LULCs selected as indicators of landscape RES and the land metrics calculated were used to describe the ecologically significant characteristics of the three types of municipalities with different degrees of coupling between supply and demand, identified across the rural-urban gradient. The analysis was performed using a mean comparison test that allowed us to characterize a qualitative category (groups of municipalities) by quantitative variables (LULCs and land metrics, respectively). For this, we used a Fisher F-test (k < 2) to determine the statistical significance of the quantitative variables in the municipality types. The more the mean of a variable in a municipality type is significantly different from the mean of that variable in the whole group of municipalities, the stronger the link between the characterizing quantitative variable and the qualitative category [58]. Those groups of municipalities characterized by a high range of statistically significant variables have greater possibilities to enjoy certain land uses and recreational services.

Steps 5. Analysis of the Relationship between Landscape Protection and the Supply and Demand for RES
The three sets of municipalities with different coupling degree detected through the rural-urban gradient by means of the product vector elements (Figure 4), comprise inside their boundaries a dissimilar representation of the PA network established in the study area. This network is composed by several conservation categories and land-use restrictions, according to national or European legislation [35]. In these municipality types, the existence of a possible relationship between the protection of the landscape and its capacity to supply the RES demanded by the NBT visitors (degree of coupling between recreational supply and demand), was tested using Fisher's F test (k < 2).

Step 6. Effectiveness of Protected Areas in Providing RES and Satisfying Visitors' Demand
The same test was used to detect, at landscape scale, the potential of PAs to provide a significant coupling between supply and demand for RES in the whole rural-urban gradient studied. To this end, we designed an inside/out approach, considering the PA boundaries and identifying the municipalities with a land area >25% inside the PA network and the neighboring ones, not included within its protection limits [1,2,59]. The coupling vector calculated (product vector; Figure 3c) allowed us to analyze the significant differences in the values of supply-demand correspondence inside and outside the PA network.

Landscape Potential to Supply Recreational Services
The quantification of the selected RES proxies as a percentage of occupied land area along the rural-urban gradient, making them spatially explicit, allowed us to characterize the potential of the landscape to provide recreational services. The data matrix collecting these proxies (Figure 3d) served as the basis for the subsequent identification and characterization, through questionnaires, of the RES of greatest interest to NBT visitors, as well as to analyse and map the relationship between the recreational supply and demand of the landscape (Figure 4a).

Visitor Preferences. Assessment of the Demand for Recreational Landscape Services
The identification of visitor preferences was based on an analysis of the frequency of responses to a set of survey questions, from which we were able to quantify the RES of the landscape with greater weight as potential tourism attractors (Figure 4b). Positive visitor responses to each of the questions posed about the recreational services provided by the landscape often exceeded 50% (10 of the 16 questions). The most valued RES were those related to aquatic systems (riparian forests, rivers, wetlands, ponds and reservoirs), with more than 80% positive responses (Table 3). By contrast, the least valued services were those linked to agricultural systems, which in no case reached 30% positive valuations. RES associated with pasture systems (grasslands and dehesas) and Mediterranean tree formations and shrublands occupied an intermediate position in the assessments, standing in a range of variation of positive responses between 42-75%. The vector of response frequencies thus obtained allowed us to establish the demand for each RES by visitors (Table 3).

Measuring the Coupling between Supply and Demand for Recreational Ecosystem Services
The spatially explicit correspondence between RES supply and demand along the rural-urban gradient was calculated by multiplying the RES supply matrix by the demand vector (Figure 4a,b). The result of this matrix operation was a vector product that allowed us to establish numerically the adjustment of the relationship between supply and demand for RES in the set of municipalities of the gradient (Figure 4c). This spatial coupling can be mapped (Figure 4d) and allows the generation of integrated RES supply-demand maps at different administrative and management scales that, in the case of this study, can vary from municipal to supra-municipal or regional scale, depending on the objectives of the land planning and management. Figure 5 represents the spatial projection of the obtained product vector, divided into three categories of degree of coupling of RES by means of the natural break classification method [53]. The resulting types of municipalities express a clear decreasing coupling gradient, in a North-South direction, linked to the rural-urban gradient, being the municipalities located to the north of the study area those that express a greater degree of correspondence between supply and demand for recreational services. Appendix A shows the municipalities belonging to each of the three types constituted according to the degree of coupling between the supply and demand.

Landscape Structure, Land Protection and Supply-Demand Coupling for Recreational Ecosystem Services
We calculated the effect of the landscape structure on the coupling between recreational supply and demand considering the composition (from LULCs) and configuration (from land metrics) of the landscape of the three types of municipalities with different degrees of coupling. Fisher's F test highlighted a significant relationship between the variation of the landscape spatial patterns and the RES coupling along the rural-urban gradient (Table 4a).
The results obtained show that the landscape of the municipalities with a lesser degree of coupling (Type 1), located to the southwest of the study region ( Figure 5), is characterized by agricultural lands, mainly rainfed crops, with a very homogeneous spatial configuration. In contrast, the group of municipalities with a higher RES coupling (Type 3), located to the north of the rural-urban gradient ( Figure 5) presents significant values of richness and fragmentation of land uses and landscape heterogeneity (PRD, SPLIT and SHDI respectively; Table 4a), which results in the spatial disconnection between patches, as indicated by the Euclidean nearest neighbour distance (ENN ; Table 4a). This allows to quantify the isolation of the patch (see Table 2). Statistically significant LULCs explain that Type 3 corresponds to a cultural landscape characterized by Mediterranean forest and traditional silvopastoral land uses. Type 2, with a medium degree of coupling between supply and demand for RES, is located in the centre and southeast of the study gradient ( Figure 5) and presents a transitional landscape, with elements typical of both agricultural lands and Mediterranean forest systems. Landscape metrics also indicate a spatial configuration with intermediate values in relation to the patterns expressed by the types of municipalities 1 and 3. Thus, the analysis of the spatial patterns of the landscape along the rural-urban gradient has allowed us to distinguish the link between the decrease in the heterogeneity of the landscape and the degree of coupling between supply and demand for RES (i.e., the less heterogeneous the landscape, the lesser degree of coupling between the supply of RES and its demand by visitors). Table 4. Characterization of the types of municipalities with different degree of coupling between recreational supply and demand, identified along the rural-urban gradient studied. Relationship between the degree of RES coupling of the municipal types and: (a) Landscape structure: (a1) LULCs, considered as proxies of RES; (a2) landscape patterns, calculated from landscape metric indices; (b) Protected area network, with different management categories and protection degree. Data are expressed as mean values per municipality types. Statistically significant values (Fisher F-test; p-value ≤ 0.05) are indicated in bold. Likewise, the analysis of the presence of PAs in the territory corresponding to the different types of municipalities detected (number of PAs and statistical significance), showed a variation gradient of land protection in North-South direction along the study rural-urban gradient (Table 4b), linked to the decrease of the degree of coupling between landscape supply and demand for RES. Thus, the group of municipalities type 3 (maximum value of coupling) is the one that has a statistically significant greater number of PAs of different categories (Red Natura 2000, Biosphere Reserve, Regional Park and National Park). The types of municipalities 2 and 3 are characterized by comprising within their limits a Regional Park and Special Protection Areas for Birds, respectively.

Degree of Supply-Demand Coupling of RES
In summary, through the rural-urban gradient studied we have observed the relationship between the degree of coupling of RES supply and demand with the heterogeneity of the landscape and its degree of protection. Thus, the territories with the greatest capacity for supplying the RES demanded by visitors (north of the gradient) are spatially more heterogeneous and have a greater landscape protection by means of different management strategies ( Figure 6). The results of the inside/out analysis designed to test (using Fisher's F test) show the effectiveness of the PA network to provide a supply of RES capable of satisfying the demand of NBT visitors. The mean value of the coupling between supply and demand for recreational services inside the protected area network is statistically significant higher than outside(150.26 and 112.04, respectivelly; p-value ≤ 0.05).

Discussion
Land use dynamics and management decisions are influenced by human demand for the services provided by ecosystems [60]. In this sense, changes in social demands have led to new trajectories of change in land use; such is the case of land preservation with the aim of conserving nature or providing recreational spaces [61,62]. Thus, multifunctional rural landscapes and peri-urban areas provide clear evidence of changes in land use driven by the demand for multiple ES [34]. In particular, new social demands, mainly of an increasingly urban population, have emerged for highly valued cultural services, such as outdoor recreation, NBT, protection of cultural heritage and different products that reflect the cultural identity of a region [17][18][19]62]. Despite this, the importance of the CES is not recognized in land use and planning schemes [63] and environmental policies do not seem to favor the participation of the population in the management of natural resources [64].
Generally, ES assessment has been performed considering the potential supply of ecosystems from a spatially explicit perspective [33]. However, few studies that assess the demand for ES (i.e., considering the amount of services required or desired by society) have modeled visitor preferences comprehensively and consistently and have produced information with spatial expression [12,65]. Furthermore, even fewer studies have analysed both supply and demand for ES from an integrated perspective [66,67]. Especially, CES are the least quantified and mapped services [68,69]. The reason for these disparities is probably related to the lack of well-established, clear and reproducible scientific methods and reference systems [60,69]. This study has among its main objectives to improve the existing methodology by making it spatially explicit.
The main contribution of the developed method is that it provides an easily replicable tool to quantify and map both the supply and demand of RES and the degree of coupling or correspondence between the set of recreational demands and the potential of the landscape to satisfy them. The method used identifies and quantifies the RES from a set of proxies based on LULCs. We selected these proxy indicators, with spatial reference, for their significant representation in the landscape and for being easily and quickly perceived by visitors to the study area [39,70]. The supply of many ES depends on their spatial context [71]. Several authors have studied the effects of the composition and configuration of the landscape on the supply of individual and multiple ES [72][73][74][75]. In this paper, we characterize the spatial configuration of the landscape using landscape metrics, which have been considered in previous analyses in the study area, as effective indicators to explore the causes and ecological meanings of landscape patterns and their consequences on ES fluxes [34,35]. The elaboration of a matrix of georeferenced data from RES has made it possible to identify and quantify the potential of the territory to provide these ecosystem services and carry out activities based on outdoor and nature tourism. The developed method facilitates the detection of possible tourist attractors on a local and regional scale.
Likewise, the demand for RES by visitors, with different perceptions and preferences, is a fundamental aspect of sustainable tourism and must be taken into account in land planning and management strategies. The survey method used was useful to determine the preferences of the visitors and their demand for RES. It is a standardized, fast, efficient and replicable technique that, considering the acceptance or rejection of landscape elements, facilitates the quantification and cartographic expression of the recreational demand [23]. The structure of the questionnaires, with a reduced number of specific questions and responses to predetermined alternatives, avoided the problems that can be caused by other approaches in which respondents must answer open-ended questions according to their own ideas. A significant number of visitors (367 people from a total of 400 contacts made) responded reliably to the survey, partly due to the structure of the questionnaire, with adequate and clear questions, easy to understand [76].
The matrix operation procedure has allowed us to quantify the correspondence or coupling between the recreational demand of visitors, based on their preferences, and the potential of the landscape to supply the required services. This spatial relationship between supply and demand for RES can be established with different degrees of adjustment [23,39]. In this work, the coupling of services has been established in a range of three categories (high, medium and low coupling), but the product vector resulting from the matrix-vector multiplication allows classifying the degree of coupling between supply and demand for RES in as many classes as observations (row vectors) compose the RES supply matrix (i.e., in this study there could be 36 coupling categories, as many as municipalities comprise the rural-urban-studied gradient). This spatial coupling can be mapped and allows the generation of maps at different administrative and management scales, depending on the objectives of land planning.
The three types of municipalities identified represent a coupling gradient of recreational supply and demand, superimposed on the rural-urban gradient studied, in the N-S direction ( Figure 5). This coupling gradient is linked to the ecological characteristics of the territory and the gradual transformation of the rural landscape, with a social-ecological dynamic that generates systems in the process of rural-urban transition [34]. The characterization of these municipalities considering proxy indicators of the landscape structure and the presence of PAs within their administrative borders, has allowed to detect the gradual variation both in the spatial heterogeneity of the territory and in its level of protection ( Figure 6). The landscape metrics used have proven to be valid and appropriate for landscape planning and landscape monitoring [77]. Thus, the municipalities located in the northern mountainous area of the rural-urban gradient (Sierra de Guadarrama and piedmont), mainly with silvopastoral land uses and the highest coupling values between supply and demand of RES in the study area, present a heterogeneous landscape with a high degree of protection. On the contrary, the agricultural municipalities placed at the southern end of the rural-urban gradient have the lowest values of these landscape indicators. The identified social-ecological pattern is strongly linked to the intensity of use of the territory, whose transformation and degradation threatens the sustainability of the ES provided.
Regarding the potential of the PA network to supply RES, the inside/out approach performed shows that the PAs studied are "hotspots of outdoor recreational services", with a greater capacity to supply the services demanded by the NBT visitors than the rest of the territory. The highest coupling values between supply and demand for RES found inside the limits of the PA network support the definition of protected area by the International Union for Conservation of Nature (IUCN) as "'A clearly defined geographical space, recognized, dedicated and managed, through legal or other effective means, to achieve the long-term conservation of nature with associated ecosystem services and cultural values" [78]. In fact, PAs are considered to offer relevant opportunities to supply cultural services, such as recreation, derived mainly from their high biodiversity and degree of naturalness [79,80]. Both of them are undoubtedly important drivers of park visits.
Nevertheless, the social-ecological representativeness of PAs and the effectiveness of their management schemes are currently highly controversial. This is mainly because PA guidelines frequently inhibit some human activities, promoting land abandonment and loss of rurality, with a negative impact on the well-being of local populations, who are vulnerable to the establishment of PAs [1,2,59,81,82]. In recent decades, it is increasingly recognized that PAs must not only focus their conservation efforts on protecting wildlife and landscape naturalness, but also maintaining and supporting the livelihoods of local people by means of more effective policies addressed to protect traditional culture, because most of the PAs of the world show some degree of human use or 'culturalness' [29,83,84].
The social-ecological approach of NBT recognizes both the need to favour the quality of life of rural people and the conservation of the resources of a valuable cultural landscape, with a great potential to expand the supply of ES [28]. Recreation in protected areas is a way of engaging people with landscape conservation measures and supporting conservation goals, beyond the appreciation of biodiversity and naturalness. However, achieving a broader social commitment with conservation of PAs requires inclusive and proactive political measures that consider a wider variety of beneficiaries and the implication of different social groups through participatory planning procedures [85,86], especially local communities and their recognized traditional ecological knowledge and cultural values [29,87]. It is noteworthy that the spatial location and the boundaries of these areas do not always coincide with the perceived landscape by local people or visitors [88]. In this regard, the designed model could be applied in an inclusive context considering both the preferences of visitors to a given territory but also the opinions and attitudes of the local population towards tourism. Thus, local population must be involved in the design and management of PAs and interested in their conservation, not only for the economic benefits of NBT, but also to achieve an environmentally sustainable development that contributes to human well-being.

Conclusions
This paper provides a social-ecological framework to face outdoor recreation and NBT as a relevant representation of CES and provides an easily replicable methodological approach that favors an integrated land management, using ecological and social resources.
The method followed has allowed quantifying and mapping, in a spatially explicit way, the spatial relationships between recreational opportunities and public perceptions and preferences, determining the degree of coupling between supply and demand for recreational services. This social-ecological interaction can be expressed at multiple scales, based on different spatial reference units. In this work the municipalities of the rural-urban gradient were used as analysis units, but the methodological approach can be expressed at different administrative and supra-municipal management scales, depending on the objectives of land planning. The analysis performed has allowed to classify the rural-urban gradient in types of social-ecological systems with different degree of coupling between supply and demand for RES and dissimilar landscape patterns as well as to detect a close relationship between degree of coupling, landscape heterogeneity and landscape protection.
The inside/out approach applied to PAs shows that they are effective hotspots of RES, with a greater capacity to supply recreational services than the rest of the territory. This proves that management guidelines of PAs and conservation-policy processes are more related to naturalness and wildlife conservation than to culturalness and the integration of local people and their cultural practices into conservation schemes. Based on these results, we emphasize the need to align nature conservation policies, tourism planning and traditional knowledge and practices, which can lead to and facilitate sustainable development.
The methodological scheme has proven to be a useful tool for strategic and sustainable land planning across rural-urban gradients and for the participatory design of PAs, generating conclusive and implementable results that are necessary in landscape management.