The various socio-economic and environmental challenges which modern society is facing nowadays, such as climate change, biodiversity loss and overexploitation of natural resources, call for the need of different systems and tools to monitor and inform the policy-makers and businesses about the development ways towards sustainability (Mooney et al. 2004, IPCC 2014 (Pachauri and Meyer 2014), IPBES 2018 (Scholes et al. 2018), IPBES 2019 (Diaz et al. 2019)). These systems and tools (Banerjee et al. 2016), embedded in academic research, offer policy-makers robust analytical frameworks, facilitating evidence-based decision-making and fostering a deeper and common understanding of the complex interplay between economic, social and environmental dimensions in the pursuit of sustainability.

There are numerous examples of such systems. In the context of climate actions, the enhanced transparency framework (ETF) under the United Nations Framework Convention on Climate Change (UNFCCC) is a reporting and review system for climate data, including GHG inventories, to track progress and provide information for policy development. The natural capital accounting represents another example of an important additional tool for providing information on sustainable development (Guerry et al. 2015). The natural capital is defined as the stock of natural resources, which includes rocks, soil, air, water and living organisms (Pearce 1988), that generate a flow of ecosystem services (Costanza et al. 2014). The framework, providing a systematic way to measure and report on stocks and flows of natural capital, is the Ecosystem accounting under the System of the Environment and Economic Accounts (SEEA-EA). It represents a valuable tool, developed by a long-lasting science-policy interface (Ekins et al. 2003, Edens and Hein 2013, Guerry et al. 2015, Costanza et al. 2017, Ruijs et al. 2018) aiming to integrate the economic and environmental data together with social data into a single, coherent framework for holistic decision-making (Edens et al. 2022). The implementation of the SEEA-EA requires diverse knowledge and expertise, such as economics, environmental science and data analysis. Accurate and reliable data, along with advanced modelling capacity, are fundamental to the effective implementation of the SEEA-EA (Hein et al. 2020).

The Ecosystem accounting includes spatial modelling of ecosystems towards organising biophysical data, measuring ecosystem services, tracking changes in ecosystem assets and linking this information to economic and other human activity through five separate accounts (Hein et al. 2015, La Notte et al. 2019, Hein et al. 2020). Like this, it presents a coherent and comprehensive view of ecosystems. The SEEA-EA framework also addresses the so called 'thematic accounting', a method that arranges data based on specific policy-relevant environmental themes like biodiversity, climate change, oceans and urban areas (Lange et al. 2022). For example, carbon accounting is part of a dedicated thematic account on climate change, but at the same time, the analysis of the carbon stock and carbon sequestration is an important part of the assessment of the ecosystem condition and services.

Soil organic carbon (SOC) stock is a critical component of the global carbon cycle, making it essential not only for carbon accounting within the SEEA-EA natural capital accounting framework, but also under the UNFCCC, where accurate calculations on the emissions due to SOC stock changes should be reported.

Estimating soil carbon content and stocks requires information on soil properties, which are measured through sampling and laboratory analyses. This is a labour and cost consuming activity (Singh et al. 2012). Soil parameter measurements on a national scale are usually part of a dedicated soil inventory or monitoring programme. In Bulgaria, the National Soil Monitoring Programme, established in 2004, is the main source of information on the relevant soil characteristics in the country. The system operates on a 16 km x 16 km grid and re-measures every 5 years. Data coming from the soil monitoring programme are used to derive country-specific emission factors to report emissions and removals from the soil carbon pool due to land-use changes as part of the national GHG inventory in line with the UNFCCC regulations. However, the current climate change policy introduces reporting-based targets, which require improvements in methodologies for reporting to achieve a higher level of accuracy. This poses a challenge for GHG reporting authorities as the application of greater precision in emissions estimates is linked to the availability of high-resolution measurement data and modelling approaches. This requires a combination of various high quality data sources in order for the calculation to reflect differences in parameters and management across different areas.

There are similar challenges when implementing ecosystem accounting at different scale – national, regional or local. As the ecosystem accounts take a spatial approach, the analysis and the assessments are presented using maps that bring together geographical, environmental, ecological and economic information in one place. If the accounting is performed at regional or national level, it is usually not sufficient to rely on individual point observations to spatially represent the soil condition and services, given the large vertical and horizontal variability of soil characteristics. Thus, large set of data as well as modelling work is required to better represent the spatial variability (Le Noë et al. 2023).

By facing these challenges, the international scientific community developed suitable tools at global level, such as SoilGrids, GlobalSoilMap.net, FAO SOIL PORTAL etc. as global sources for soil information. They are developed, based on different types of data and modern technologies. One of the main advantages of these systems is that they provide spatial distributions of soil properties across the globe. This makes the datasets suitable for analysis and spatial representation of the ecosystem condition and services under SEEA-EA. In addition, these tools have the potential to contribute to better estimates of the country-specific emission factors in the GHG Inventories under UNFCCC, especially when the available dataset does not provide enough measurements to average out the sampling error, when grouping soil measurements into climate regions and soil types.

SoilGrids is a system for digital soil mapping, based on global compilation of soil profile data and environmental layers. Digital mapping methods are constantly evolving and improving; particularly popular nowadays are machine-learning models and algorithms, which, according to Khaledian and Miller (2020)and Poggio et al. (2021), provide better results than most multiple linear regression methods. However, the application of these tools in regional and national assessments requires that model data be validated with independent information, despite their acceptable levels of accuracy (Poggio et al. 2021).

In this regards, the aim of the current study is to compare the modelled data from SoilGrids with regional soil information by presenting a methodological approach for verification of the suitability of SoilGrids data for the purpose of ecosystem reporting under the natural capital accounts or the GHG inventory under the UNFCCC regulations.

Case study area

Methodological approach

The suitability of the SoilGrids data for ecosystem reporting has been verified against the independent set of data. The correspondence of the predicted values from SoilGrids has been compared to randomly distributed point observations from the dataset. The comparison has been done at plot level and then at sample average. The parameters which have been compared are pH, bulk density, coarse fraction, soil organic carbon (SOC) content, nitrogen content and soil organic carbon (SOC) stock.

The SoilGrids data are downloaded in a raster image format at a spatial resolution of 250 m for each of the studied indicators (pH, org. C, N, bulk density, coarse fraction, C stock) and by layers – 0-5, 5-15, 15-30 cm. SoilGrids data on soil C stock are presented on the web-based platform only for the 0-30 cm soil layer. By using a mask layer of the study area, the analysis is guided to the range of assessment. Raster calculator was used to derive data for each indicator for the 0-30 cm depth for subsequent analysis and processing (Fig. 1).

Figure 1.

Conceptual scheme of the study.

To verify the accuracy of the model data, we used an independent set of data from a few sources - ICP Forest plots in the region, published data from scientific study and our own observations (Suppl. material 1). The soil measurements from these sources differ in comprehensiveness and layers surveyed (Table 1). This makes validation by layer depth difficult. In order to reduce the error resulting from filling in missing information (e.g. for bulk density at depth and SOC content for 20-30 cm), the comparison of model results was performed for 0-30 cm depth and not by layers. The only exception was the bulk density, which was validated for 0-5 cm, as the model data correspond exactly to the depth analysed in the field using a ring-type metal probe.

Table 1.

Download as

CSV

XLSX

Soil parameters by source of information and soil depth.

Depth	0-5 cm	5-15 cm	5-20 cm	0-10 cm	10-20 cm	15-30 cm	20-30 cm	0-30 cm
SoilGrids (raster data, 250 m resolution)	pH, CF, BD,SOC content, SOC stock, N content	pH, CF, BD,SOC content, SOC stock, N content				pH, CF, BD,SOC content, SOC stock, N content		SOC stock (derived)
ICP Forest (soil profile data)	BD			pH, CF, SOC content, N content	pH, CF, SOC content, N content		pH, CF, SOC content, N content	SOC stock calculated according to IPCC GPG LULUCF (2003)
Eli Pavlova-Traykova (2018) (soil profile data)	pH, CF, BD,SOC content, N content		pH, CF,SOC content, N content					SOC stock calculated according to IPCC GPG LULUCF (2003)
Own observations (soil profile data)	BD			pH, CF, SOC content, N content	pH, CF, SOC content, N content		pH, CF, SOC content, N content	SOC stock calculated according to IPCC GPG LULUCF (2003)

We analysed the agreement between model data and field observations by parameters both at the sample plot level and as sample average. The averages of the two sets have been tested by two-tailed paired samples t-test and Wilcoxon signed-rank test.

To facilitate the comparative analysis further, ancillary materials, such as climatic data, soil types and terrain models, were incorporated into the estimation process to reflect stratification needs and improve the efficiency of the analysis. The forest vegetation zoning is used as a proxy for the combined effect of these factors. According to the Bulgarian classification on forest vegetation zones, there are three main vegetation belts in the region of Rila Mountains - The Lower Plains and Foothills oak woodland (0-700 m a.s.l.), The Mid-montane beech and conifer forests (701-2000 m a.s.l.) and the zone of High Mountain area (> 2000 m a.s.l.).

Data sources

Laboratory analysis

Calculating SOC stocks (0-30 cm)

The analysis of the main physico-chemical properties of the soils in Rila Mountains, based on the data of SoilGrids, shows that the bulk density of soils ranges from 0.93 to 1.45 g/cm3. The mean value was determined to be 1.18 g/cm3 (± 0.1). Soil bulk density variation in the study area exhibits a clear trend of decreasing with increasing elevation (Fig. 2).

Figure 2.

Soil bulk density in the Rila Mountains for 0-30 cm depth. Source: SoilGrids. The Barplot is a histogram with x: bulk density (cg/cm³) and y: frequency.

In terms of coarse fraction, the soils in Rila Mountains are characterised as grainy (Fig. 3). The coarse fraction varies between 4.5% and 30%. Most often it is around 14.5% (± 4.4%). The higher the altitude, the greater is the coarse fraction. The value above 2000 m a.s.l. varies between 20 and 30%.

Figure 3.

Coarse fraction content of soils in the Rila Mountains for 0-30 cm depth. Source: SoilGrids. The Barplot is a histogram with x: coarse fraction (cm³/dm³) and y: frequency.

The pH (for 0-30 cm) is mostly slightly acidic to neutral (pH 5.5-6.5 and above 6.5). The mean pH (H₂O) was determined to be 6.03 (±0.4) according to SoilGrids data. From (Fig. 4), it can be seen that there is a clear pattern of increasing soil acidity with increasing elevation.

Figure 4.

pH values of soils in Rila Mountains for 0-30 cm depth. Source: SoilGrids. The Barplot is a histogram with x: pH and y: frequency.

The content of the soil organic carbon (for 0-30 cm) in the soils of Rila Mountains varies between 5.1 and 96.7 g/kg (Fig. 5). Values between 20 and 60 g/kg predominate, with the mean content determined to be 38.5 g/kg (± 13.35). The nitrogen content of soils ranges from 1.44 g/kg to 5.60 g/kg (Fig. 6). The mean value of nitrogen (N) concentration is calculated to be 3.21 g/kg (± 0.9). In both parameters – C and N concentrations, the content of both SOC content and N is higher with an increase in altitude.

Figure 5.

Soil organic carbon (SOC) content in Rila Mountains for 0-30 cm depth. Source: SoilGrids. The Barplot is a histogram with x: SOC content (dg/kg) and y: frequency.

Figure 6.

Nitrogen content of soils in Rila Mountains for 0-30 cm depth. Source: SoilGrids. The Barplot is a histogram with x: N content (cg/kg) and y: frequency.

When comparing the information from SoilGrids with the data from field observations in the Rila Mountains, it is noticeable that the results from the field studies show a greater variation in soil parameters compared to those predicted by the model (Table 2, Fig. 7). This is to be expected taking into account that there is high variability in soil indicators even within relatively homogeneous soils and environmental factors (Wells and Case 1995). However, with respect to the mean values of the individual indicators, comparable results are reported for SOC content, coarse fraction, bulk density and satisfactory with respect to N content and pH values (H₂ O) (Table 2).

Figure 7.

Distribution of values of the studied parameters by type of data.

The comparative analysis and the statistical test performed (Student t-test; paired, two tailed, p = 0.05) showed that, in terms of bulk density data, coarse fraction, SOC content and C stock, no statistically significant differences were found between the sample means from the field data and those from the models (Table 3). The tests show that there is a significant difference between the means of the field observations and the model for pH and N.

As the assumptions for normality of the data from the field studies are not completely met and there are also outliers in some of the analysed parameters, an additional non-parametric test was performed to check the statistical significance in means - the Wilcoxon Signed Rank test (paired, two tailed, p = 0.05). The results confirmed the outcomes from the t-test (Table 4).

The SoilGrids data for the Rila Mountains represents well the main soil parameters in relation to their depth and dynamics in terms of altitudinal range. The analysed parameters take into account the known variability resulting from the complex influence of a number of climatic and biological factors (Franzluebbers et al. 2001, Sergeevna Kozun et al. 2022), which are closely correlated with changes in altitude and tree cover (Figs 2, 3, 4, 5). It is noticeable that the content of org. C and total N increase with increasing altitude (Table 5), since the accumulation and decomposition of organic matter are directly dependent on the hydrothermal regime of the soil and the vegetation and tree composition, which strongly differentiate with altitude. Soil acidity also varies with altitude from acidic to slightly acidic in the lowlands to strongly acidic in the mid-montane and high-montane zones. Bulk density decreases with increasing altitude because of the increase in soil porosity. This is due to accumulation of organic matter and formation of thick and soft forest litter, contributing to the formation of thick humus horizons. Coarse fraction increases with height, especially in the brown forest soil’s zone, which is generally characterised by high skeletal content (Koynov et al. 1998) and is significant in the alpine zone of the mountains, where there are mostly cliffs.

With regard to nitrogen content, it is striking that model and field data have significant differences in values (Fig. 7, Tables 3, 4). More than half of the reported data from field observations fall outside the model predicted minimum nitrogen content values. In SoilGrids, nitrogen content includes total nitrogen (ammonia, organic and reduced nitrogen) as determined by Kjeldahl, plus the absorbable forms of nitrogen, nitrate and nitrite. The method we applied in the analysis of the field observations does not account for nitrate and nitrite forms of nitrogen, but their content in the forest soils is very small (Donov et al. 1974) and the difference cannot be explained by a difference in their accounting alone.

SoilGrids data predict pH changes well with increasing elevation and changing vegetation cover, with significantly more acidic soils in the conifer forest zone (Fig. 4). However, the comparative analysis shows that, overall, the model predicts less acidic soils with pH values mostly between 5.7 and 6.6, while field observations report pH values between 4.5 and 5.8 (Fig. 7). In the Rila Mountains, Cambisols are the most common soil type (> 50%). Data from several studies on these soil types in the country confirm the more acid solution reaction (< 6.0) (Zhiyanski 2018, Kirova 2019, Malinova et al. 2020, Hristov et al. 2021). The acidity of soils in the Rila Mountains is conditioned by the widespread occurrence of carbonate-free silicate rocks - granite-gneisses, on the one hand and the predominance of coniferous forest tree species on the other hand, which lead to soil acidification (Binkley and Sollins 1990, Zhiyanski 2018).

We suggest that the reason for the discrepancies between model and field data in terms of pH and nitrogen content is likely due to model imperfections associated with the input data. The data in SoilGrids are modelled on the basis of information from nearly 240,000 soil profiles from around the world, dominated by data from agricultural land. In addition, it should also be taken into account that more than 60% of these data were collected between 1960 and 2000, 34% have an unknown sampling date and only 16% of soil profiles were surveyed between 2001 and 2020 (Poggio et al. 2021). Batjes et al. (2020) note that the sampling period should be taken into account when mapping using environmental covariates, especially when this refers to highly variable features that may be affected by land-use changes or natural disturbances. However, Poggio et al. (2021) considered that the spatial variation of dynamic features is greater than their temporal variation and this was neglected in the creation of SoilGrids.

The analyses and comparisons with field data show that SoilGrids provides reliable values for soil organic carbon content, bulk density and coarse fraction, which are important parameters for determining the amount of carbon stored in the soil. The SOC stock is a dynamic characteristic that is influenced by land use (Zhiyanski et al. 2016, Menichetti et al. 2017) and silvicultural practices (Ameray et al. 2021) in addition to the soil potential itself. The analysis of the average soil organic carbon (SOC) stock by soil type across the Rila Mountains corresponds to the observed SOC stock in forest soils at national scale reported by Stoeva and Kirova (2021). According to SoilGrids data, intersected with the soil map of the area, the mean SOC stock in Cambisols in Rila Mountains is 57.92 tC/ha, while the average SOC stock in Luvisols is 53.00 tC/ha. Stoeva and Kirova (2021) reported mean SOC stock values in interval 51-56 tC/ha in Cambisols and 43-53 tC/ha for Luvisols. Both studies show that the SOC stocks are characterised with large variability. It is highest in brown and dark-coloured forest soils, as well as in upland-meadow soils, which are typical for mountainous regions in Bulgaria (Hristov and Filcheva 2017). The reason for the large dispersion in the values of soil parameters (including SOC) is explained by the diverse conditions under which soil formation processes take place. Brown and dark-coloured forest soils in Rila Mountains develop in a moderately wet and cool climate under the influence of woody vegetation - mostly moist beech and conifer forests. The relief is hilly to mountainous and, in places in the Rila Mountains, it is highly dissected, which determines the development of erosion processes. Depending on the altitude, the different types of forest vegetation and the geological composition, certain differences in the characteristics of individual brown soils are noticeable, which is also clearly visible from the large dispersion in the SOC content in brown and dark forest soils (Fig. 7). The mountain-meadow soils, on the other hand, are developed under harsh climatic conditions, contributing to the formation of a peaty horizon, water retention and slowing down of the mineralisation of organic matter. This leads to humus accumulation, delayed chemical weathering, strong leaching and acidification.

Regarding the SOC content, the most significant influence is that of vegetation and mineralisation processes. In this respect, woody debris is the main source of organic matter, nitrogen and ash elements in the soil beneath forest ecosystems. Climatic factors in turn determine the processes of mineralisation and transport of organic matter and ash elements in depth. The cycling of substances under beech, spruce and alpine shrub formations is different, which has implications for the properties and characteristics of soils in the mountains and can be clearly seen by looking at the soil indicators and their variation with height, which is followed by the variation in vegetation and climatic conditions. When comparing soil organic carbon stock (org. C) data from SoilGrids with different forest vegetation zones (according to Bulgarian classification), there is a clear difference in average values and spread across various altitudes. The lowest forest zone has low C stock, while the high mountain zone has higher C stock (Fig. 8). There is no great difference in average soil C stock between the sub-zones at High Mountain area (above 2001 m a.s.l.), but there is a noticeable difference within the average SOC stock values between the sub-zones of the other two forest vegetation zones – the one of the Lower Plains and Foothills oak woodland (0-700 m a.s.l.) and the Mid-montane beech and conifer forests (701-2000 m a.s.l.). The difference in the mean SOC values of the Oak woodland zones is explained mainly by the influence of the land use and its legacy effect. Most of the agricultural land is present at the lowest altitude 0-500 m a.s.l., where the lowest average soil C stock is estimated. Regarding the difference in means of the sub-zones of the Mid-montane beech and conifer forests, the variation is explained by diversity of the environmental factors within this zone, such as terrain, climate and vegetation. Cambisols are the main soil type distributed within this zone. Based on the SOC analysis by soil type, Cambisols have the highest dispersion in SOC values. The average SOC value within the sub-zones differentiates substantially. These findings are consistent with the observations of Kyuleva (1979), who proposed the division of brown forest soils into three subgroups. The first subgroup shows transitional characteristics with the grey and cinnamon forest soils (Luvisols), the second consists of the typical brown forest soils (Dystric-Eutric Cambisols) and the third shows transitional characteristics with the dark-coloured forest soils (Umbric Cambisols).

Figure 8.

SOC stock distribution and density by vegetation zones and its altitude ranges. SOC data: SoilGrids.

SoilGrids data show that the SOC stock is also heavily influenced by the local climate, specifically temperature and moisture. This is evident when combining the SOC data from SoilGrids with information on climate. The highest soil carbon storage is observed in the cool and humid climate zone, which prevails in the Rila Mountains. This climate promotes the accumulation of organic matter in the topsoil, particularly fulvic acids (Hristov and Filcheva 2017). The humidity is also a relevant factor for SOC stocks in the investigation area. Even within similar temperature ranges, there is a significant difference in the mean soil carbon stock between cool, dry areas and cool, humid areas (Fig. 9).

Figure 9.

Distribution of the SOC stock data by climate zones (IPCC). SOC data: SoilGrids.

In warm, dry climates, soil carbon storage is low and humic acids dominate the organic matter. This is exemplified by the eroded and leached forest soils, which are typical in these climate zones and have some of the lowest SOC levels in the Rila Mountains.

The analysis confirms the importance of considering climate regions, soil types and vegetation when deriving country- or region-specific SOC stocks and soil carbon stock changes for ecosystem accounting purposes, especially when calculating greenhouse gas inventories under the UNFCCC. To achieve this, stratification should be applied to group soil measurements by climate region and soil type, ensuring that each stratum has sufficient data to minimise sampling error. When using international databases such as SoilGrids, rather than averaging all pixels or plots within a region, more precise approaches, such as averaging plots with similar climatic and soil conditions - even from neighbouring regions or countries - or using statistical methods to identify the most comparable plots, can provide more accurate results (Bellassen et al. 2023).

SoilGrids data represent a valuable resource for assessing key indicators of soil health and SOC stocks, providing critical information for different ecosystem reporting domains, such as the Ecosystem Accounts (SEEA-EA) and GHG inventories. One significant advantage of SoilGrids is its accessibility and spatially explicit data. This makes it useful for the representation of the soil parameters in the ecosystem accounts, which are usually represented by maps.

SoilGrids aims to provide soil information on a global scale. However, its application at regional or local scales requires careful consideration. Independent verification of SoilGrids data with field observation from Rila Mountains (Bulgaria) demonstrates that SoilGrids data can be effectively utilised for regional analysis and reporting in Bulgaria within the natural capital accounting framework or for other reporting needs. However, special consideration should be given to the prediction of the N content and pH in mountainous area, as the study also demonstrated.

SoilGrids is also a vital resource for research and analysis when soil information is scarce, offering a means to fill data gaps. However, as a model output, SoilGrids possesses inherent limitations and uncertainties. Therefore, instead of relying on specific point data, it is advisable to conduct analyses at a broader scale, focusing on areas with similar environmental conditions (e.g. climate, elevation, vegetation, soil type) to those lacking data. This approach helps mitigate the uncertainties associated with the model while still providing meaningful insights.

The research that obtained these results was supported by the LTER-BG infrastructure (Agreement No. DО1-320/30.11.2023), purchased under the National Roadmap for Research Infrastructure, financially coordinated by the Ministry of Education and Science of the Republic of Bulgaria.

LS: Conceptualisation, Methodology, Analysis, Visualisation, Writing original draft. MZ: Writing, review and editing.