The study area is comprised within the Central Balkan National Park and encompasses the mountainous regions of Stara planina. The park partially covers territories of the administrative regions Gabrovo, Lovech, Plovdiv, Sofia and Stara Zagora. The park boundaries enclose the resort of Beklemeto. The region of the Central Balkan National Park Reserves includes, partially or fully, the wildlife sanctuaries Boatin, Dzhendema, Kozya Stena, Peeshti skali, Sokolna, Stara reka, Steneto, Severen Dzhendem and Tsarichina, where 28% of the land is recognised as areas with a strict protection regime.

In landscape-ecological terms, the area is characterised by rich forest biodiversity dominated by broad-leaved forests of Querqus cerris, Querqus frainetto, Fagus sylvatica and mixed forest patches with Pinus sylvestris (Cudlín et al. 2017). About 70% of the ecosystems are natural. There are 17 forest habitats designated in the list of special protected areas of the EU Habitats Directive.

In the studied area of the Central Balkan National Park, forests and natural vegetation cover 51.35% of the territory, where 37.55% is covered with broad-leaved forests; coniferous forest represent 0.69%; mixed forests occupy 6.80%; and natural grasslands – 3.24%. The agricultural activities are mostly developed in the northernmost and lower parts and nearly half of the territory of the study area (47.0% or 556124.00 ha) is used for cultivation of various crops and permanent pastures (Fig. 2). Still intensely anthropogenised and built-up, the natural vegetation is preserved and the broad-leaved forests dominate the land cover within the region.

The main concerns of the study area are related to naturally and anthropogenic-driven hazards triggered by climate change, by natural disasters (avalanches, straight-line winds, fires etc.) and by land intensification. Forestry and agricultural practices confront the demand for habitat restoration and conservation to sustain the overall health of the environment. Considering the importance of the region in terms of Green Infrastructure (GI) and biodiversity, the scientifically-based evidence and the integration of the ES concept to monitoring information is expected to raise the awareness of both the scientific and the decision-making communities and to support the ongoing environmental projects for park management.

Every ecosystem holds a specific potential to deliver ES defined by its functions. The socially driven demands for these services are the factors that transfer an ecosystem function into a service (Fig. 1). The integration of ecological indicators that provide information about the environmental status, conditions and interconnections between the landscape components provide opportunities for assessment of the present ecological condition and trends (Grêt-Regamey et al. 2015, Kandziora et al. 2013, Maes et al. 2016, Müller and Burkhard 2012, Tammi et al. 2017, Zhiyanski et al. 2009, Zhiyanski et al. 2008). For the realisation of the objectives of the present research, the authors apply the ecosystem services concept that allows the assessment and mapping of the dynamics in soil properties in forest ecosystems.

The empirical analysis, related to the potential services provision, highlights the ecological status of the forests in the Central Balkan region. The biophysical indicators correspond to the definition of services and are derived from biophysical data provided by the ICP Forests monitoring data for the case study area. The indicators define the “final” benefits that could be provided by the ecosystems within a particular land cover unit. According to CICES V 4-3, these indicators also serve as a tool for quantitative assessment and mapping of certain ecosystem services such as Decomposition and fixing processes and Global climate regulation by reduction of greenhouse gas concentrations (Table 1).

The indicators proposed for the identification of the ES from CICES - pH, thickness of the organic layer, soil type, concentrations of carbon and nitrogen, are further developed and supplemented with more details in the MAES framework (Maes et al. 2016). However, the available records from the surveys carried out at ICP Forests sites provide comprehensive information for assessing the above-mentioned services. Therefore, it was decided to combine the suggested indicators with the defined indicators concerning green infrastructure in the framework of the project MetEcosMap.

The ICP Forest monitoring network detects and keeps record of laboratory analysis and information on the soil profiles for each site (ICP level I) including data for pH (CaCl₂), Total Organic Carbon (g/kg), Total Nitrogen (g/kg), mean dry bulk density of the soil, heavy metals etc. Later, these indicators were used to calculate the Soil Organic Carbon (SOC) stock and the Soil Organic Matter (SOM) as indicators for the provision of regulating ecosystem services.

Soils have a crucial role in the ecosystem functioning. The soil type is an important landscape indicator and can be used to characterise the forest ecosystems. Therefore, it was important to identify the different soil types for better characterisation of the mapping units. Spatial information for the soil types was obtained from the JICA Project database and was adapted to FAO-WRB.

The ecosystem integrity and the landscape peculiarities depend on the cumulative interaction of all natural and social components described in the human-environmental systems. The CLC 2012 data obtained from the Copernicus land monitoring inventory was used to delineate polygons as the smallest-scale spatial units to identify and map the ES. Given the land cover dataset, relevant up-to-date information can be used as a reference base for further investigations.

The data processing and the elaboration of maps was realised under the ESRI ArcGIS platform. The application of the ecosystem services concept is relevant to depict the spatial and temporal dynamics of the soil properties in the forest ecosystem functioning. As the study concerns especially forest lands (incl. forest stands and lands with natural vegetation), ecosystem types such as urban areas, water bodies and bare rocks have been excluded from the selection.

At the next stage, the spatial overlay of the two sets of data (CLC and soil) was made to generate spatial units for analysis, assessment and mapping of the ecosystem services (Equation 1).

Equation 1: CLC12/(Soil Type) = UNIT ID

Once the polygons have been identified (e.g. polygon with unit ID #053 is characterised by mixed broad-leaved forests on Planosols), they were connected to the plots of the existing network – ICP Forest Level I. The specific objective of this step was to be able to observe the changes and the potential of the forest ecosystems to provide ecosystem services with regards to the dynamics of the soil properties. To demonstrate the highest potential to supply the particular ecosystem services, two periods were considered that can detect the dynamics in forest ecosystems. The available records in the ICP Forest monitoring network date back to 1992/93 and the last ones are from 2015. Therefore, the period from 1992/93 until 2015 has been delineated and used as a historic reference to track the changes to the present.

The mapping approach provides information about the spatial distribution and allocation of the defined services. One of the main advantages of this method is the opportunity to represent and visualize data collected from field surveys and empirical estimates in GIS. The approach for mapping and assessment in forest ecosystems, developed in the national methodology, was applied to identify the appropriate indicators and to rank their values with scores in order to describe the capacity of the ecosystem to supply selected regulating ES. Based on information and data collected from the ICP Forests monitoring network, biophysical data have been estimated and have been associated with each delineated spatial unit. The main goal was to gather quantitative data and to map territories representative for the region, where the anthropogenic influence is relatively low. The selection of quantitative indicators and their integration in the mapping process uncovers the actual ecosystem service supply of forests. At a later stage, the values have been transferred into GIS and have been aggregated to a 5-point Likert scale (1-5) considering the assumption that the interval between the values is equal and that each interval contains values.

Statistical distribution was applied to the entire ICP data sample. Firstly, the Total Organic Carbon (TOC) and the Soil Organic Carbon (SOC) values were classified for each soil layer and assessed with a corresponding score (1-5). Later, when assessing the ecosystem services, the estimates were considered for the entire soil profile (0-20 cm), which resulted in reclassification of the given ranks for the separate indicators.

The data, derived from the analysis of the organic carbon in soil for the period from 1992/1993 to 2015, illustrate the trend in the parameter values (Fig. 3). Corresponding to the intervals for the ES assessment, the entire sample is ranked between 1 and 5 (soil layers 0-10 cm and 0-20 cm), where the mean equals to 32.32, the median is 27.75 and the coefficient of variation is calculated to 0.80. The majority of the values is clustered in the interval 0-40, representing an equivalent score 1 – very low. In addition, a certain decrease in the estimated TOC values is observed. For observation year 2015, few monitoring forest sites can be ranked with medium score - 3 and there is none with a higher score for the entire study area of the Central Balkans.

To demonstrate the highest potential to supply the particular ecosystem services, two periods are considered that can detect the dynamics of the soil properties in selected forest areas in the region. Data for different indicators collected from ICP Forests Level I plots were analysed for the period from 1992/93 (first observation) to 2015. The carbon stored in the topsoil layer (soil organic carbon) is identified by two variables. The first is the mean dry bulk density of the fine earth in kg/m3. The second is the TOC contained in the soil organic matter of each soil type in layers from 0-10 cm and 10–20 cm. The values were used to rank this indicator and to associate with the relevant score 1-5. The SOC (t/ha) is given by the following equation:

Equation 2: Soil organic carbon (SOC)=C×H×Q×K

Where:

C – carbon content (%), H – thickness of the soil layer (cm), Q – bulk density (g/km3), K – (100 – Ck)/(100) (%) correction coefficient considering skeletal fraction Ck.

The SOC represents 0.58% of the soil organic matter stock. Therefore a factor of 1.72 (Van Bemelen factor) was used to convert from SOC to SOM [%]. The estimates resulting from the carbon content in the topsoil were further used as indicators for the ecosystem service Decomposition and fixing processes. An increase in SOM stock and in the C concentrations is thus recognised as leading to greater biological diversity in the soil and to having regulating functions over the biological control of plant diseases and pests. The Level I data principally describe the topsoil layer to a depth of 20 cm. In general, the decrease in SOC content is proportional to the increase in depth for mean SOC content and is estimated below 20% (Hiederer 2009). However, the lack of measurements recorded in the database for the soil parameters implies incomplete soil profile data and certain limitations.

Statistical distribution to define the score was also applied to the SOC estimates. The mean, median and coefficient of variation, based on the data obtained for the entire time of observation, were calculated. For the period 1992/93-2015, the median is 38.25 and the mean equals to 44.06, giving a coefficient of variation of 0.67, which characterises low variance in the distribution of SOC data for the period of 20+ years.

Later the data were distributed into equal intervals and this method was used to classify the values into a certain number of categories with ranges equal to 40 in 5 classes with attribute values from 0-140. The datasets for the ICP Level I observation plots were distributed in 5 categories (ranks) from 1-5 (Fig. 4).

The decrease in the overall SOC concentration for the 2015 dataset indicates higher SOC accumulation over time, which equals the C sequestration rates. The climatic conditions (rainfall and temperature), soil characteristics (clay concentration), soil erosion and degradation and land use practices (arable land, grassland, agricultural farming) are some of the potential variables related to the analysis of the SOC stock. Interestingly, for observation year 2015, most of the values are clustered in the intervals 10-40 and 40-80 and, for observation year 1992/93, they vary between 6.08 and 38.42 t/ha.

The SOM stock, in combination with the GI area, is also one of the representative indicators for the ES Global climate regulation by reduction of greenhouse gas concentrations. Globally, as containing over two-thirds of the soil carbon, forest ecosystems are considered as significant carbon sinks that play an important role in the carbon cycling and in climate regulation. The GI index is derived from calculations on the basis of amount of GI per capita. In the present research, for the Central Balkan region, the GI index is estimated on the average value of the available indices for the adjacent municipalities (Table 2). The index is calculated in hectares per capita (haPC) and can be successfully applied as it refers to the area of natural and semi-natural areas aggregated from the LUISA land use map.

The functional relationship between the biotic components in forest ecosystems can also be characterised with the processes frequently affected by changes in biodiversity including climate regulation and carbon sequestration. To perform a comprehensive analysis of the ecological dynamics, however, complex factors corresponding to the status of forest ecosystems, ecological cycles and terrain aspect as well as man-made influences, should be considered.

The estimates demonstrate that, for polygons with natural grasslands (code 321, CLC 12) on Planosols, there are no changes detected in the last 20+ years (Fig. 5). Interestingly, the tendency that can be observed is that the state of the forest ecosystems is deteriorating as the quantities of soil organic matter are decreasing and the area of the territories, where the period 1992/93 is outlined with medium supply capacity (score 3), in 2015 is characterised with low supply of that service and score 2. These changes are tracked within areas occupied by broad-leaved forests.

Estimations for TOC and for GI per capita indicate significant variations in broad-leaved and mixed forests on Planosols and Cambisols. For a period of 20+ years, the score corresponding to the capacity for supply of this service within broad-leaved forest has dropped from low to very low (Fig. 6).

According to data for observations from 2015, only natural grasslands on Planosols have better capacity to supply that service. For the rest of the territory, it can be assumed that the dynamics in the soil properties in selected forest ecosystems have had a negative effect.

The Supply matrix on Table 3 shows that, for the Central Balkan study area, Planosols with natural grasslands and Cambisols with coniferous forests have the highest potential to supply ES from the group of Regulation & Maintenance (CICES V 4-3). The range is defined in scale 1-5, where 1 is very low supply capacity and score 5 is very high supply capacity.