The MAES methodological framework provides typology for ecosystems, a set of indicators for the assessment of ecosystem condition and mapping of ES (Directorate-General for Environment (European Commission) 2013). The core elements of the framework are: 1) mapping of ecosystems; 2) assessment of ecosystem condition; 3) mapping and assessment of ecosystem services. The main steps which cover the core activities of an operational framework are presented in the conceptual scheme developed through the ESMERALDA*1 project (Burkhard et al. 2018a). Furthermore, within the project, a multi-tiered approach for ES mapping and assessment was developed. The approach considers different methods (biophysical, socio-cultural and economic) at different levels of detail and complexity and can be applied according to specific needs, data and resource availability (Burkhard et al. 2018b). However, every ecosystem assessment has to be relevant to a certain theme and address a broad range of questions pertaining to decision-making processes that occur at different levels of decision-making and across different actors of society (Burkhard et al. 2018a). The mapping and assessment of ES provided by the NH for the needs of recreation and tourism is a theme that could contribute to finding solutions to the problems related to the conflicts between conservation and the use of the NH. A conceptual framework for mapping and assessment of ES provided by the NH in Bulgaria for recreation and tourism was developed through several efforts within the framework of the Centre of Excellence “Heritage BG”*2. It is based on the assumption that the generation of NH for the needs of tourism can be presented as the linkages between the natural systems and tourism in the form of ES potential, flow and demand (Nedkov et al. 2021b). A set of indicators for mapping and assessment at a national level are proposed and the methods for their quantification are arranged following the multi-tiered approach. They allow the production of ES maps for the priority ES which can be used for planning purposes in sustainable tourism. ES maps quantify and visualise where and to what extent ecosystems contribute to human well-being (Burkhard and Maes 2017). To represent ES provided by NH in a spatial context, it is necessary to define where ES are generated and what is the potential of the ecosystems. To map the overall potential of an area, it is necessary to integrate the whole range of ES. The matrix approach enables normalising all ES values in a uniform score from 0 to 5 which makes possible the integration of the resulting GIS layers.

Study area and initial data

The multi-tiered approach, developed in this work, is designed for application at the national level. Therefore, the whole area of Bulgaria is selected as a case study. 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 in terms of species present on the territory, it hosts 26% of all European species, 70% of the protected bird species under the EU Birds Directive and 40% of the conservation habitats types (under Annex I, Council of the European Union 1992). Both biodiversity and geodiversity, as elements of NH, are major sources of recreation and tourism in the country (Ihtimanski et al. 2020). The country is divided into nine tourist regions (Fig. 1). Each of them has a city recognised as a centre of the region and priority tourism branches which make the specialisation of the region. For instance, the Rhodopes Region is specialised in mountain and religious tourism, while the Valley of the Roses is in health and cultural tourism (Nikolova et al. 2021b).

Figure 1.  

Case study area. Tourist regions: D - Dunav (The Danube); DR - Dolina na Rozite (Valley of the Roses); RP - Rila and Pirin; R - Rodopi (The Rhodopes); SC - Severno Chernomorie (North Black Sea coast); S - Sofia; SP - Stara Planina (The Balkan); T - Trakia (Thrace); YC - Yuzhno Chernomorie (South Black Sea coast).

MAES implementation needs spatially-explicit datasets to address the key drivers, pressures and their different gradients and variations in space and time. Each ES is assessed by specific indicators which have to be supported with appropriate spatial data available at the national level corresponding to the whole territory of the country. The activities under MAES in Bulgaria led to the development of several datasets, but their use at the national level at this stage is hampered by two main problems. Firstly, the data for the nine ecosystem types are in separate datasets which do not fit topologically correctly if they are merged in a single GIS layer. Secondly, the mapping does not cover the Natura 2000 areas which is a significant gap that makes these data inappropriate for national scale mapping. Only for some services, which are assessed using municipality-based initial information, there are appropriate data that can be applied at the national level. For instance, the quantification of education and science service is based on a number of papers calculated per municipality (Assenova et al. 2018).

The lack of full coverage at the national scale data of some ES can be overcome using models and modelling approaches. For instance, the ESTIMAP model for recreation uses easily available data on land cover, protected areas, water bodies, transport network and topography (Ihtimanski et al. 2020). The various data sources used for ES modelling at tier 3 are presented in Table 1. Nevertheless, there are still some services that could not be quantified using available data at the national level. Expert-based assessment is the only possible option for these services. The matrix approach at tier 1 necessitates appropriate spatial units to be selected for the assessment. In this study, we use the map of ecosystem subtypes, based on CORINE land-cover (CLC) data correlated to the MAES ecosystems classification at the third level for Bulgaria (Hristova and Stoycheva 2021). It contains 27 ecosystem subtypes which can be delineated after the correlation with CLC classes.

The multi-tiered approach for mapping and assessment of ES provided by the NH at the national level is based on the MAES framework (Directorate-General for Environment (European Commission) 2013) and on the conceptual scheme for its implementation (Burkhard et al. 2018a). It consists of four main elements (Fig. 2) which are described in the following subsections.

Figure 2.  

Conceptual scheme of the multi-tiered approach. The four boxes in the scheme correspond to the four subsections below.

The application of the multi-tiered approach enabled us to develop a GIS database containing layers for each of the 15 priority services (Suppl. material 2), as well as integrated layers about the overall potential and the potential of the main ES (provisioning, regulating and cultural). The GIS layers were used to prepare maps of ES provided by the NH in Bulgaria at the national scale (Fig. 3A-C). The main results present the overall potential of ecosystems in the country (Fig. 3D) which is the most important output directed to various practitioners in the recreation and tourism activities. The map shows that the areas with very high potential are almost evenly distributed across the country. They are scattered in polygons of various sizes. However, several clusters with a concentration of polygons with very high potential areas could be outlined. The largest one is located in the south-western part of the country within the high mountain areas of Rila, Pirin and Western Rhodopes. The second one covers areas in Central Stara Planina (Balkan Mountain) and Sredna Gora Mountain. There are also three other clusters in the Stara Planina located in its western and eastern parts. Two clusters are formed in the low mountain and hilly areas of Eastern Rhodopes and Strandzha. The last one is located in the lowland-hilly area of the Eastern Danube plain. The areas with very high potential cover 9578 km² which is about 9% of the country.

Figure 3.  

Potential of the natural heritage to provide ES. A - provisioning, B - regulating, C - cultural, D - overall. 0 - no potential; 1 - very low potential; 2 - low potential; 3 - moderate potential; 4 - high potential; 5 - very high potential.

The areas with high potential cover about 24% of the country (Table 4) and their distribution shows a similar pattern as the very high potential. They form compact areas in the Western Rhodopes and Central Stara Planina Mountain and more or less scattered areas in the other mountains. The areas with moderate and low potential are predominant covering together almost half the country. The areas with very low potential are located predominantly in the lowland areas in the northern and south-eastern parts of the country comprising about 15% of its area. In general, the results show that almost the whole country has some kind of ES potential provided by the NH and only 4% are assessed as with no potential.

The maps of provisioning, regulating and cultural services visualise quite different patterns of ES potential throughout the country. The overall potential of the provisioning services (Fig. 3A) is quite low compared to the other two groups of services. The areas with higher potential for the provisioning services are located in the northern and south-eastern parts of the country where the topography is predominantly flat and the agricultural ecosystems are widespread. Regulating services have just the opposite pattern with higher potential in the southern and central mountainous areas (Fig. 3B). Cultural services have higher overall potential with slightly higher values in the mountains (Fig. 3C).

The ES assessment of the NH for recreation and tourism enabled also the estimation of the ES potential per tourist region. We recalculated the overall ES scores for each tourist region estimating an average ES score as well as the distribution of the 0-5 scores. The average scores show quite similar results for all regions with figures ranging from 2.45 to 3.08 (Fig. 4). Two predominantly mountainous regions (Rila and Pirin and the Rhodopes) have the highest scores which exceed 3.0. In the third place is the Valley of the Roses region which comprises both mountainous and lowland areas in the central part of the country. The two Black Sea coast regions are in the “middle of the table” compared to the other regions. The Danube and Trace regions that occupy the areas with lower potential have also lower average scores. However, the region with the lowest average score is The Balkan which covers mainly mountainous areas.

Figure 4.  

Potential of the tourist regions to provide ES. A - Mean potential of the tourist regions to provide ES, B - Distribution of the potential scores within the tourist regions by area. Tourist regions: D - Dunav (The Danube); DR - Dolina na Rozite (Valley of the Roses); RP - Rila and Pirin; R - Rodopi (The Rhodopes); SC - Severno Chernomorie (North Black Sea coast); S - Sofia; SP - Stara Planina (The Balkan); T - Trakia (Thrace); YC - Yuzhno Chernomorie (South Black Sea coast).

In contrast to the relatively uniform average scores, the distributions of the 0-5 scores amongst the regions show quite different patterns. Each region has specific distribution and only the first two (Rila and Pirin and the Rhodopes) show a similar pattern in the distribution diagram with high and very high potential covering more than half of the area, moderate and low potential covering the rest, while the areas with 0 and 1 score have limited extent. The Danube Region has predominantly low and very low potential as they cover more than 75% of the area. The North Black Sea coast has a similar pattern with a slightly higher share of the area with moderate potential. This is in contrast with the South Black Sea coast which has a significantly higher share of the areas with high and very high potential than the North Black Sea coast.

The results for the ES potential were obtained using various methods at three different tiers. From the methodological point of view, it is important to compare the results at different tiers. There are only two ES assessed by methods at tier 2 which is not enough for appropriate conclusions. Therefore, the analyses were made only for tier 1 and tier 3. The mapping results (in the form of GIS layers) were re-arranged into two groups corresponding to these tiers. The layers were processed to recalculate the ES potential derived from the method at different tiers. At tier 1, the scores for the nine ES from the expert assessment were recalculated to estimate mean values for each ES. Then, the mean values were normalised to the 0 to 5 assessment scale. At tier 3, there were seven layers produced by the different modelling methods. They were processed using the same procedure which was performed for the integrated layer of the overall ES potential. Thus, we had two resulting layers representing the results about ES potential, calculated using the methods at tiers 1 and 2. These scores could not be treated as another way to define the potential of the NH to provide ES. They are just for analysing the results at different tiers and to obtain data for discussion about their advantages and disadvantages from the methodological point of view. This enabled us to generate maps of the ES potential derived from methods at tier 1 and tier 3, as well as the differences between them (Fig. 5).

Figure 5.  

ES potential and the differences at tier 1 and tier 3. Tourist regions: D - Dunav (The Danube); DR - Dolina na Rozite (Valley of the Roses); RP - Rila and Pirin; R - Rodopi (The Rhodopes); SC - Severno Chernomorie (North Black Sea coast); S - Sofia; SP - Stara Planina (The Balkan); T - Trakia (Thrace); YC - Yuzhno Chernomorie (South Black Sea coast).

The two maps of the ES potential show a similar pattern which correlates relatively well with the overall ES map presented in Fig. 3. However, the differences in the actual scores are pronounced between the two maps. The areas of lower ES potential in the map of tier 1 are assessed mainly as 2 (low potential), while in the map of tier 3, their scores are mainly 1 (very low potential). The areas with no potential (score 0) are better represented on the map of tier 3, while on the tier 1 map, they are limited to smaller patches. The areas with very high potential are larger on the tier 1 map, while on the tier 3 map, they are smaller and more fragmented. Furthermore, areas of very high potential can be found also in areas where the tier 1 map indicates lower potential. The distribution of the areas with different ES potential is given in Table 5. It shows that the results from the two tiers coincide well only for the areas of moderate ES. The areas of very low and high potential are significantly higher in the tier 3 results, while the areas of low and very high potential are higher in the tier 1 results.

The comparison between the results obtained by methods at tier 1 and tier 3 were analysed by overlay between the two layers. First, the scores at tier 3 were recalculated to negative values. Then, an overlay procedure by a simple adding operation between the two layers was applied. Thus, in the areas where the scores are equal, the resulting value would be 0, in the areas where the tier 1 score exceeds the tier 2 score, the result will be a positive value between 1 and 5 depending on the excess value, in the areas with a higher score for tier 3, the result would be a negative value with the same gradient. The result of this procedure was a new layer presenting the differences in the scores between tier 1 and tier 3 (Fig. 5). The results show that, in 91% of the area, the scores are equal or differ by one unit (Fig. 6). The areas with excess higher than two units are negligible (below 1%). In general, the scores at tier 1 are higher as the areas where they exceed by one unit cover 42% percent of the studied area. The coincidence between the two scores is found for 36% of the area. The areas with an excess of tier 3 scores are located mainly in the mountains, while the tier 1 excess is mainly in the lowlands (Fig. 5).

Figure 6.  

Comparison between tier 1 and tier 3. A - Spatial distribution of a negative image between both layers where T3 exceed T1 and vice versa, B - % of the whole area.

Data availability and accuracy of the resulting ES maps are amongst the most important issues in the application of the tiered approach (Grêt-Regamey et al. 2015). The accuracy of the maps strongly depends on spatial resolution and the quality of both initial and intermediate data sources. Intermediate data sources in our study are the data layers generated as the result of different spatial analyses during the indicators’ quantification stage of the study (see Fig. 2). The analysis of the spatial resolution of the data sources at different tiers would bring important information about the accuracy of the resulting maps and the uncertainty analyses of the ES scores. The main data source for assessment at tier 1 is the CLC dataset which was processed to generate the ecosystem subtypes layer. This layer contains the spatial units used in the expert assessment and it was used to analyse the data at tier 1. The statistical information at tier 2 is available at the municipality level; therefore, the municipality’s outlines should be considered as spatial units for the mapping. The modelling approaches applied at tier 3 combine different data sources and, as a result of the modelling procedures, the output layers have different data characteristics. To compare the quality of the data at the different tiers, we combined the layers from the assessment of the different ES. The comparison of the data quality is made using four characteristics of the polygons from the resulting vector layers (Table 6). The spatial resolution of the resulting data at tier 3 is the highest, while at tier 2, it is much lower. The resolution at tier 1 is in the middle of the others.

The results of the ES layers generated at tier 3 significant differences in the spatial resolution (Table 7). The highest spatial resolution has the resulting data from the LCZ model applied for the assessment of the local climate regulation (ES VIII). The lowest resolution has the results from the assessment of the aesthetic value (ES IX) by the InVEST model. The results from ESTIMAP about recreation are comparable to the spatial proxy methods applied for regulation of natural hazards (ES VI) and maintaining population and habitats (ES VII).

The assessment of multiple ES at the national level is a challenging task because it necessitates a variety of data that should be available for the whole country and the application of various methods that requires a large team of experts with different expertise. This is possible only for large and well-funded projects that are not easy to be achieved. Even the ES matrix (which is easy and not resource-intensive) is more often used at the local and regional than at the national level (Campagne et al. 2020). The multi-tiered approach presented in this study ensured the mapping and assessment of 15 ES at a national scale and the integration of the results for the needs of one specific activity (tourism). The dataset which was developed during this study enables easy and fast generation of national scale maps for various purposes that can ensure further analyses to support the sustainable use of the NH for various tourism activities in the country. Particular examples are the studies on specific tourist activities, such as outdoor tourism (Ihtimanski et al. 2020), ski tourism (Silvestriev et al. 2021), forest therapy (Dodev et al. 2021) and speleological tourism (Nikolova et al. 2021b).

The results demonstrate that the NH of Bulgaria is a valuable resource that ensures the generation of various ES which are important for the development of tourism activities in the country. The areas with very high potential can be found throughout the country which proves the hypothesis behind tourism regionalisation which covers the whole country and distinguishes the regions depending on their specialisation. The clusters of very high potential correspond to six out of nine tourist regions. The Rhodopes Region contains two clusters and has also one of the highest overall ES potential scores. Stara Planina is the other region with two clusters, but its overall score is lower. The reason behind this difference could be explained by the more compact mountainous character of the Rhodopes Region and the high forest cover. Both mountain relief and forests cover stand out as the main factors for the high ES potential. Thus, the Stara Planina Region contains also some lowland areas with a higher anthropogenic impact which reduces the overall score of the region. Although the increase in the elevation tends to refer to an increase in the ES potential, the highest areas in Rila and Pirin are not assessed with the highest potential. In this case, the lack of forest in the alpine and subalpine areas is the factor for the decrease in the overall potential. This could be defined as one of the limitations of the approach that needs to be studied in more detail in the future. The application of some kind of a rapid assessment approach that exploits available datasets and triggers more detailed and disciplined specific studies on ecosystem condition indicators (Kokkoris et al. 2018) is needed.

In this work, we develop and apply an approach for mapping and assessment of the NH as a source of ES for recreation and tourism. It is based on the multi-tiered approach proposed by Burkhard et al. (2018b) which relies on integrating different methods at different levels of detail and complexity. The authors state that it “can be applied to specific needs, data and resources availability”. The multi-tiered approach presented in this paper builds on the above-mentioned approach by developing its application one specific need, the assessment of NH for recreation and tourism. The results of its application at a national level in Bulgaria prove its applicability and potential to solve such a complex task. The main advantage of the approach is the possibility to assess various ES with different data availability and specifics which necessitates the application of different methods. This enables the development of NH tourism assessment beyond the usually assessed cultural services (such as outdoor recreation, cultural heritage and aesthetic experiences) to the relevant regulating (such as maintenance of habitats and local climate regulation) and provisioning services (such as water supply and crop production). Thus, the assessment of the ES potential for tourism can be better related to nature conservation which is crucial for the preservation of the NH and the achievement of sustainable tourism.

The integration of the ES matrix into the approach allows for the assessment of more ES, especially at tier 1, as it helps to overcome the limitations of data availability and the lack of proper proxies for quantification (Campagne et al. 2020, Prodanova and Varadzhakova 2022). At tier 3, we apply a combination of already proven and well-documented modelling tools such as InVEST and ESTIMAP with new modelling approaches developed for this study. The application of the LCZ model for the assessment of local climate regulation has a high potential as it relies on the freely-available Copernicus dataset and a tool for LCZ delineation which is applicable to all EU countries (Demuzere et al. 2019). The spatial proxy modelling approaches developed for two ES (water for drinking purposes and maintaining of population and habitats) rely on spatial analysis techniques available in the main GIS software packages (such as ArcGIS and QGIS) and easily-available spatial data. Thus, they can be applied in various areas for studies with limited resources.

The multi-tiered approach was applied predominantly at different levels of scales. The most representative example is provided by Grêt-Regamey et al. (2015) by tier 1 at the continental level, tier 2 at the national level and tier 3 at the local level. In our study, the three tiers are applied at the national which gives the opportunity to compare the results obtained by methods at different tiers. The comparison between the results from the expert-based assessment at tier 1 and the more complex modelling methods at tier 3 shows differences that rarely exceed 1 unit on a 0 to 5 scale. Therefore, the results at these two tiers have a general agreement and the modelling results can be considered as a validation of the expert assessment. The deviations from this general agreement show that the expert-based assessment gives slightly higher scores. One possible explanation could be in some kind of exaggeration of the NH by the experts. They have been asked to evaluate the potential of the NH and this term refers to something valuable. Furthermore, the experts evaluate ecosystems defined from land-cover data. For instance, they give a very high score to the deciduous forest and this score is transferred to all such forests in the dataset, which is an inherent limit of the expert-based matrix as discussed in many papers (Jacobs and Burkhard 2017, Campagne et al. 2020). However, the modelling methods at tier 3 rely on more indicators that reveal the inherent heterogeneity of the forest and different scores according to this heterogeneity.

The spatial data resolution at the different tiers can be used as an indicator of the data quality and consequently of the accuracy of the results. The spatial resolution of the tier 2 data is quite low due to the specifics of the data which is available at the municipality level. The multi-tiered approach could be further developed by considering the specifics of the ecosystem types, especially the necessity of finer-scale mapping of urban and freshwater ecosystems. As Haase et al. (2014) mentioned, scales appropriate for urban ecosystem analysis need to be developed. On the other hand, mapping at a national scale necessitates easily-accessible open data instead of the more precise, but resource-intensive unmanned aerial vehicle data. The Copernicus data provides appropriate means and new opportunities for the assessment of urban areas (Sarafova 2021). There is also a lack of more detailed data for some ecosystem types, such as sparsely vegetated and shrublands. Recent habitat mapping studies (Grigorov et al. 2021) provide some data, but further studies are needed.