Mapping the natural heritage as a source of recreation services at national scale in Bulgaria

Natural heritage includes natural features or natural areas of outstanding universal value. At a national level, this value refers to the importance of ecosystems which can be considered as the spatial units representing the natural heritage of the particular area in terms of their values to people. Nature-based outdoor recreation represents an important service that interests millions of people and contributes to connecting them to nature, but it may also cause negative impacts in the form of pollution, erosion and habitats loss. We apply the ESTIMAP recreation model which provides a framework for a spatially-explicit assessment of local outdoor recreation and use it to identify and assess the natural heritage as a source of recreation services at a national level in Bulgaria. At the first stage of the study, we identify the natural heritage and the data sources to represent it in a spatially-explicit way. Then, we apply the module for recreation potential to assess the potential of the natural heritage to provide a recreation ecosystem service. At the third stage, the accessibility of the natural heritage is assessed in order to specify how the potential identified at the previous step can be really used. Finally, the recreation potential and accessibility are integrated into the recreation opportunity spectrum in order to develop the maps representing the ecosystem service supply provided by the natural heritage. The results are presented in form of a recreation potential map that reveals the capacity of natural heritage to provide the recreation potential, map of the accessibility of the natural heritage and map of the recreation opportunity spectrum representing the combination ‡ ‡ § © Ihtimanski I et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. between the first two maps. The maps will be used for the development of an innovative geospatial platform designed to facilitate the access of the Bulgarian natural heritage to the European common knowledge and innovation markets. The results on the accessibility and recreation opportunity spectrum contribute to the development of the model in areas which were not covered by previous applications at the EU scale.


Introduction
According to the World Heritage Convention, Natural Heritage (NH) includes natural features consisting of physical formations, geological features and physiographical formations, natural sites or precisely delineated natural areas of outstanding universal value from the point of view of science, conservation or natural beauty (UNESCO 1972). Transposed at a national level, the outstanding universal value is converted as "natural significance", which refers to the importance of ecosystems, biodiversity and geodiversity for their existence value, as well as their scientific, social, aesthetic and life-support value (Harrison and O'Donnel 2010). Ecosystems incorporate biotic and abiotic elements (i.e. biodiversity and geodiversity) and can be considered as the spatial units which can represent the NH of the particular area in terms of their values to people. These values are usually described as the ecosystem services (ES) delivered by nature to people (MA 2005). The NH can be related mainly to cultural services, such as outdoor recreation, tourism, cultural heritage, aesthetic value, but also to some regulating services, such as maintenance of habitats and local climate regulation. Nature-based outdoor recreation, as an ES, has been studied and spatially defined using various methodological approaches. Most commonly used methods combine GIS models with participatory approaches applied in different scales -at a local level (Nahuelhual et al. 2013;Kienast et al. 2012;de Vries and Goossen 2002), at regional level (Tardieu and Tuffery 2019;Kliskey 2000;Koniak et al. 2011;Birch et al. 2014;Garnache et al. 2018;Yuxi et al. 2018;Sonter et al. 2016), at a national and supranational level (Willibald et al. 2019;Balzan and Debono 2018;Vallecillo et al. 2019).
Outdoor recreation represents an important service that interests millions of people and contributes to connecting them to nature. It includes both local and long term recreation. The former is not considered as tourism because there is no overnight stay, while the latter can be represented as touristic outdoor recreation. These activities have an important role in human well-being and health since they provide physical, aesthetic and cultural benefits and offer an opportunity to experience directly a relationship with nature (Zulian et al. 2013;Norman et al. 2010;Lankia et al. 2015). On the other hand, outdoor recreation may cause negative impacts on ecosystems through pollution, intensified erosion, harm to wildlife or habitats and biodiversity loss (Wong 2004;Tavares et al. 2012). A response to these problems could be the concept of sustainable tourism which aims to balance the environmental, economic and socio-cultural features of tourism development by maintaining environmental resources, the socio-cultural livelihoods of host communities and providing stakeholder benefits (Schloegel 2007). The assessment of ES can be an important tool in the process to bridge the conceptual gap between the ecological and social sciences, by linking the state of ecosystems with human well-being and activities (Drius et al. 2019. The sustainable use of the NH for recreational purposes in Bulgaria is the main problem addressed by this study. It is a part of a broad research project (Center for Excellence "Heritage BG"*1) that aims to develop a methodology to promote the access of the Bulgarian natural heritage to the European Digital Single Market of Knowledge and Information Services (*2). The development of a geospatial platform that will enable an inventory and evaluation of the ES related to recreation is one of the main objectives of the project and outdoor recreation is the pilot service to be evaluated. Therefore, we need a tool to identify and assess the capacity of natural heritage to provide recreational services in a spatially-explicit manner. The ESTIMAP recreation model provides a framework for a spatially-explicit assessment of local outdoor recreation (Zulian et al. 2013;) which could appropriately meet our needs. The ESTIMAP recreation model is based on "Advanced multiple-layer Look-up Tables" which assigns ES scores to land features according to their capacity to provide the service . The values of ES scores for each input are based on literature or expert input and the final value is based on cross-tabulation and spatial composition derived from the overlay of the different thematic maps. The model is developed for the territory of EU, but its application is possible also at the local scale for different spatial units , such as protected areas (Schägner et al. 2018) and urban areas (Baró et al. 2016;Cortinovis et al. 2018;Maes et al. 2019), as well as at regional supranational scale (Liquete et al. 2016). There is still no evidence on application at the national level for a particular country and this study also aims to fill this gap. Another challenge of the current work is the interpretation of the linkages between NH and recreation through the ES concept. In the study, we adapted the original model configuration to our needs, by using a vector instead of raster data, increasing the spatial resolution, introducing topography as an indicator (following the suggestion of Schägner et al. 2018) and other aspects which are discussed in more detail further in the paper.
In this respect, the main objective of this study is to utilise the ESTIMAP recreation model to identify and assess the NH as a source of recreation services at a national level in Bulgaria. In this paper we aim; • to identify NH at a national level in Bulgaria, • to assess and map the recreation potential (RP) of the NH, • to assess and map the accessibility of the NH and • to reveal the spatial aspects of the potential benefits to people by the recreation opportunity spectrum (ROS).

Material and methods
We utilised the ESTIMAP recreation model as a methodological basis of our study ( Fig. 1). At the first stage, we identify the NH and the data sources to represent it in a spatiallyexplicit way. In this procedure, we take into account the criteria for NH identification and the accessibility of the spatial data at the national level as input for the ESTIMAP recreation model. At the second stage, we apply the ESTIMAP module for RP to assess the potential of the NH to provide recreation ES. The third stage includes the assessment of the accessibility of the NH, which allows us to specify how the potential identified at the previous step can be really used for recreation purposes. At the fourth stage, following the model workflow, the RP and accessibility are integrated into the ROS in order to develop the maps representing the ES supply provided by the NH. The locally-adapted version of the ESTIMAP recreation model is a little bit different than the original version , Zulian et al. 2013. For instance, some components or subcomponents are excluded from our adaptation, such as the distance from the coast, the proximity of inland protected areas to the coast and the assessment of benefits, while other elements, such as the topography, are included (more details in the next subsections). For all the data processing, calculations and overlay operations ArcGIS Desktop 10.2 and its tools were used.

Study area and data
The approach, developed in this work, is applied at the national level in Bulgaria, therefore the mapping was implemented for the whole area of the country. Due to the diverse climatic, geological, topographic and hydrological conditions, Bulgaria is amongst the richest countries in Europe in terms of biodiversity and geodiversity. The country accounts for about 2.5% of the total EU area, but there are 26% of all European species, 70% of the protected bird species and 40% of the conservation habitats*3. Both biodiversity and geodiversity, as elements of NH, are major sources for recreation and tourism in the country. The national ecological network consists of protected areas (under the Protected Areas Act -PAA*4) and protected zones (under the EU Birds and Habitats Directives). The protected areas are grouped into six IUCN categories and consist of 1015 sites which comprise 5.3% of the country's area (Fig. 2). There are 234 protected habitat zones (30% of the area), while the protected bird zones are 119 (22.7%).
Bulgaria has a total of 10 cultural and natural sites of exceptional value for humanity in the World Heritage List (MOEW 2019*3). Only three of them are in the "NH" category -Srebarna Lake, Pirin National Park and Old-growth beech forests in the Central Balkan National Park. The latter is included in a serial site with trans-boundary significance "Ancient and Primeval Beech Forests of the Carpathians and Other Regions of Europe".
An important aspect of NH is the naturalness of the landscapes in the country. The main source of spatial data about the landscapes is the CORINE Land Cover (CLC). It has a 3levels hierarchical classification system with five classes at the first level and 44 classes at the third level. At the first level, Bulgaria has all five classes distributed as follows: agriculture areas (51.66% of the country territory); artificial surfaces (4.79%); forest and semi-natural areas (42.55%); water bodies and wetlands (0.99%) (Fig. 3). Most of the protected areas in the country fall into forest and semi-natural areas.
The ESTIMAP recreation model needs geospatial data, which are used in various GIS and tabular operations for land features identification. Different ES scores need to be assigned to the land features according to their capacity to provide the service. The ES scores given to each land feature are based on expert assessment or literature review (Zulian et al. (2013). GIS maps are used to represent the final value which is based on zonal statistics for the land features and the ES scores received. The data used for this study are obtained from open data sources. All the data and its sources are presented in Table 1. Most of the data sources are the same as in the ESTIMAP: ES mapping at a European scale , Zulian et al. 2013). Most of the data used in our research are geospatial, only the data for bathing water quality are tabular, but it is processed and transformed into a GIS point layer. All the geospatial data we use are in vector format, except several intermediate raster layers obtained during the data processing.

Identification of the NH as a source of recreation ES
The natural elements that can be recognised as a heritage for the purposes of the study are selected according to the national classifier of the sites of importance for NH. This classifier has been developed within the framework of the project "Heritage BG"*1. It is substantiated in four main categories, of which we use: world-established NH sites and components of geo-heritage. They are selected according to their national importance and data availability. The world established NH sites include UNESCO sites and natural protected areas which are available in the CDDA database (Table 1). The components of geo-heritage include various abiotic elements, such as water bodies, caves, rock formations etc. For this study, we selected the water bodies because they can be supported by data at the national level. In addition, we selected land cover and elevation data in order to define the degree of naturalness and to represent the topography.
These categories have been combined with the ESTIMAP criteria for RP to define the following four main components for identification of NH as a source of recreation ES: Table 1.
Input data for the adaptation of the model.
The DN can be described as the difference between the current and the potential natural state of a particular area to define to what extent it was transformed by human impact. In our case, the DN is a measure to identify the areas with preserved NH from those which are transformed into systems dominated by anthropogenic elements. The latter could also have some components of the NH, but they are not recognisable when the mapping is at a national level. The most appropriate source of spatial data is CLC which is available for the whole country and its classes can be easily related to specific categories representing the DN.
The presence of NP is used as additional information to emphasise the importance of the natural areas which have already been recognised as areas with special value. This value has been legally anchored by the specific restrictive regime of each protected area. Thus, the protection status of the different protected areas can be used to categorise them according to their natural value and, respectively, to identify their contribution to the NH. We use the DN and NP as the basis to distinguish the areas with preserved NH from those without.
Water is one of the main components of NH, especially in relation to recreation (Kakoyannis and Stankey 2002). The existence of water bodies affects the biodiversity and the attraction of the area which increases the value of the NH. It also depends on the water quality as it is one of the main factors for outdoor recreation. We use two parameters to represent the water component: coastal areas and bathing water monitoring points.
The topography is also an important component of recreation as the areas with mountainous relief are usually more attractive for people (Schägner et al. 2018). Thus the mountain areas can be defined with higher NH value and the topography is used as an additional component in the assessment of the RP (see Fig. 1). In this study, we defined the elevation above 1000 m as a limit to outline the areas with a higher value due to their mountainous relief.

Assessment of the RP
The assessment of the RP of the NH is based on the four main components described in the previous subchapter which have been used to derive indicators representing the capacity of the areas to provide outdoor recreation (Fig. 1).

Degree of naturalness component:
For assessing the DN, we use CLC as a geospatial input. We use the hemeroby index as an indicator for the naturalness of the land cover classes. Firstly, the original scale of the hemeroby index is inverted in order to fit the scale used in our approach. Thus, the classes with higher anthropogenic impact receive a lower score, while those with better preserved natural elements receive a higher score. For instance, the continuous urban fabric (class 111) has the highest score (7) according to the original scheme and it is converted into the lowest (0). The scores of the CLC categories and their correspondence to the hemeroby indices are given in Table 2. A higher score is assigned to the water bodies compared to the actual value of the hemeroby index because of the assumption that they provide some special attraction for recreation activities.

Natural protection component:
For this component, we use "Nationally designated areas (CDDA)" as a geospatial input. The scores, based of the assessment presented in Zulian et al. (2013), were used as reference values which were transformed into integer format instead of floating points in order to facilitate the further data processing (

Water component:
As described above, the water component is assessed using two parameters, coastal areas and bathing water. In order to generate a geospatial input for the coastal areas, we use a CLC layer with the coastline of Bulgaria. The coastline is transformed from line to polygon (1 km buffer). For the second parameter, we use the bathing water quality status data provided by EEA in .xlsx format for the year 2018. The initial data are georeferenced and a point layer is created. Then the point layer is transformed into polygons (400 m buffer) and the scores are assigned according to Table 4. Finally, the coastline polygon and the bathing water polygon are intersected to generate WC layer.

Topography:
The layers representing the components described above (DN, NP and WC) are united in one feature class. The scores for each component are summed for any CLC class, designated area or water body. For the areas above 1000 m, the summed score is increased by 1 in order to reflect the additional component topography. Then the final score is classified according to the scheme given in Table 5.

Accessibility of the NH
At the third step of the analysis, the accessibility of the NH is addressed in order to assess how the recreation can be delivered to people. The accessibility is based on the distance of the NH to settlements and roads. The Accessibility map is generated from vector feature classes -Urban areas and Road network in the Republic of Bulgaria. The Euclidean distance is computed for the two layers. The rasters with computed distance are combined using the ArcGIS raster calculator to create the Accessibility layer. The Accessibility layer is converted to a vector map and, in the cases where a polygon falls in two zones, the higher one is taken. The polygons are assigned a score according to Table 6. Scores used for bathing water quality Table 5.
Final scores for RP based on DN + NP + WC + Topography.

Recreation Opportunity Spectrum
The ROS represents how people can benefit from the opportunities provided by nature for recreational activities if they are able to reach them (Zulian et al. 2013). In our case, it is used to map the different degrees of the recreation provided by the NH in the country, based on its accessibility for the potential end-users. The ROS has been generated by intersecting between the RP and the accessibility map. The RP has been classified into five categories and the accessibility map has been classified into three categories following Zulian et al. (2014) and . After the intersection, the final score for ROS is assigned according to the assessment of RP and the accessibility parameters as given in Table 7.  Final scores for accessibility based on roads and settlements distances Table 7.

ROS
ROS according to Accessibility and RP.

Recreation potential
The RP map (Fig. 4) represents the potential of Bulgarian NH to provide the ES recreation.
The results show that the most common areas of the country are those assessed as medium potential. They comprise about 45% of the Bulgarian territory. These are mostly areas overlapping with land categories, such as forest and semi-natural areas, agriculture areas, wetlands and water bodies. Some artificial surfaces have also been identified with medium RP and they are mainly sport and leisure facilities, discontinuous urban fabric and green urban areas.
Three percent of the country's territory is covered by areas with no potential for recreation. Most of these areas are artificial surfaces excluding green urban areas and sport and leisure facilities, which have higher potential. The areas identified with a score of 1 -Low potential are mostly agricultural areas and some forest and semi-natural areas covering around 52% from the country.
The areas with high and very high RP are less than one percent, they are mostly forest and semi-natural areas, agricultural areas and wetlands. Some sports and leisure facilities also have high and very high RP. Most of the areas with the highest RP are located on the

RP of the NH in Bulgaria.
Black Sea coast and in the high mountain areas, but they are too small and could not be easily recognised on the national scale map (Fig. 4). The location and extent of some very attractive recreation sites with high and very high RP on the Black Sea coast of Bulgaria are presented in Fig. 5.

Accessibility of the NH
The map of the accessibility of the NH (Fig. 6) shows that the settlements in Bulgaria and the roads connecting them form a dense network which covers the entire country. About 80% of the area is located within 5 km distance from a settlement or a road. Therefore, these areas are classified as "near" in the accessibility map. The second zone defined as "proximal" comprises about 17% of the country's area. The distance to roads and settlements in this zone is between 5 and 10 km. The third zone, defined as "far" has a limited extent and covers only 3% of the area. It is located mainly in the upper parts of the mountains. The largest extent of this zone is in the south-western part of the country along Rila and Pirin mountains, as well as the western parts of the Rhodopes. The central part of the Stara Planina chain represents the second largest area with remote NH. The only lowland area which falls into this class is located on the western part of the Danube plain in the area with extensive landslides which disrupt the roads network.  ROS of the NH in Bulgaria.

Recreation opportunity spectrum
The map of the ROS in Bulgaria (Fig. 7) shows that the predominant part of the country (47%) is in areas with low RP and easy access. They are located in the whole country, but predominantly in the elevation belt below 1000 m. About 30% of the country's area is with medium RP which is easily accessible. These are predominantly low mountain areas located in the southern and western parts of the country. The third place is for the combination of medium potential -accessible, which comprises 12% of the country located in higher mountains Rila, Pirin, the Rhodopes and Stara Planina. These three combinations altogether cover almost 90% of the country's area while the other six have relatively limited extent with only 10% altogether.
The distribution of the different ROS combination is given in Table 8. As is mentioned above, a limited area (less than 1%) is located in zones with high and very high potential. The places where these two zones match with easy accessibility are located predominantly at the Black Sea coast and few of them in the high mountains (Fig. 8). The places with high RP combined with easily accessible and accessible areas can be found as small patches on the map, located mainly in the Rhodopes and the Stara Planina mountains, as well as at the Black Sea coast. Table 8.

ROS
Recreation Opportunity Spectrum.

Discussion
In this study, we apply the ESTIMAP recreation model to identify, assess and map the RP of the NH at the national level in Bulgaria. The original model configuration has been adapted for the needs of the study by increasing the spatial resolution using local vector data instead of raster data and adding one new component for the assessment of the RP. This component is based on the assumption that mountain regions have a special attraction for recreation and generate higher values per visit. Such an assumption has been also discussed by Schägner et al. 2018. We also give a higher score for all the freshwater reservoirs and lakes because they provide a special attraction for nature-based outdoor recreation such as fishing, boating, swimming, rafting etc.
The ESTIMAP recreation model was applied for the whole EU, but the results for Bulgaria and Romania were included only in the map of the RP (Zulian et al. 2013). The comparison with our map of the RP shows quite a similar pattern with higher potential in the mountainous areas and lower in the lowlands. Our results on the accessibility and ROS could be considered as an original contribution to the development of the model to all European countries. The comparison of ROS distribution between the maps presented in Zulian et al. (2013) and  and our map presented in Fig. 7 shows quite similar patterns in the distribution of the low potential with easily accessible and Location of the areas with higher ROS.
accessible areas, but the high potential combination is quite different. The results from our study present a limited extent in Bulgaria, while in the rest of Europe, they have far larger extent. It is a matter of further analyses to find the reasons for such differences, but we can suppose some differences in the data sources and the scale as possible factors. The data about water bodies and the elevation are with a higher resolution which results in higher differentiation of the areas with higher RP becoming smaller. The use of elevation as RP criteria decreases the score of the lowland which have almost no areas of high RP. The use of more detailed roads network data leads to limited low accessibility areas. Another reason is the lower number and area of water bodies in Bulgaria compared to western and northern Europe.
The approach of the ESTIMAP model is intuitive and relatively easy for application at various scales. Most of the required data are available for free which facilitates its applicability. The approach is flexible to the addition of other indicators and data that can improve the quality of the results. On the other hand, data manipulation is sometimes redundant and time-consuming. The scores of the RP are somehow subjective as they are based on expert assessment, therefore further improvements of the approach could be in the form of incorporation of more empirical data for this process. Some CLC classes (322 -Moors and heathland; 323 -Sclerophyllous vegetation; 331 -Beaches, dunes, sands; 332 -Bare rocks; 334 -Burnt areas; 421 -Salt marshes; 422 -Salines; 521 -Coastal lagoons; 523 -Sea and ocean) are not included in the hemeroby index assessment given by Szilassi et al. (2017). The incorporation of these classes in the analyses would benefit further studies on outdoor recreation, especially in the areas where such classes have higher significance in the landscape pattern.
One limitation of the approach is that there are some other factors for the RP of the NH, such as climate condition, terrain morphology, fragmentation of landscapes etc. which are not incorporated in the current version of the model. There are also NH sites in areas not identified as such, for instance in urban areas. These limitations are admissible for studies at a national level due to the lack of appropriate data or the degree of their importance, but for future studies at a local level, more comprehensive adaptation of the ESTIMAP model is necessary.
The results of this study will be used as a basis to develop the ES component of a geospatial platform for access to the Bulgarian NH. Such a platform can be used for geospatial analysis of multiple threats to and from tourism (Drius et al. 2019) and the outdoor recreation assessment is a good example to test the possibilities for such analysis. The importance of outdoor recreation and its links to the natural-based tourism has increased significantly over the past few decades (Pröbstl and Haider 2013) and it is in close relation to the concept of sustainable tourism and green economy (Pan et al. 2018). The outdoor recreation assessment integrated into the platform will enable the end-users to apply analyses on both economic and nature protection aspects in order to plan sustainable tourism activities.

Conclusions
With the adaptation of the ESTIMAP model, we have managed to identify the spatial extent of the NH in Bulgaria at a national scale, to assess its RP and accessibility and derive ROS. These results could be considered as a basis for assessment of the capacity of Bulgarian NH to provide recreational ES. The maps of the RP, the accessibility and the ROS will be used for the development of an innovative geospatial platform designed to facilitate the access of the Bulgarian NH to the European common knowledge and innovation markets. The outputs of the platform will be used in the spatial planning and implementation of the recreational industries and will support the capitalisation of the future products in relation to various business entities. The results of the outdoor recreation could be used in analyses of the RP and its accessibility for the needs of the district tourism management plans. The ROS results will be integrated into the geospatial platform in order to develop a tool for place-based analyses which can be used by the tourism companies in their investment activities.
More indicators, such as climate conditions, landscape pattern and topography, as well as more detailed and precise data about the current indicators, should be added for further development of the model in Bulgaria in order to achieve more reliable and representative results. An indicator of the riparian vegetation should be added to reveal this important recreational aspect. The accessibility part should be improved with actual data on the population density because the available data are from 2011. The road network data should be developed with railways, small roads and hiking trails. The implementation of the model at a local level assumes upgrading the current data with more detailed NH objects, such as the open geospatial database OpenStreetMap etc. (as, for instance, Cortinovis et al. 2018, Cortinovis andGeneletti 2018). More extensive mapping studies on a large scale would reveal the overall picture of the potential of the NH for recreation at a national level in a socio-environmental context.