Towards soil bulk density maps with quantified uncertainty in Mediterranean soils

Soil Bulk Density (BD) is a key physical property integrating soil texture, structure, and soil organic carbon (SOC), and it exerts strong control on porosity, water and solute transport, root penetration, and gas exchange. Because BD links soil mass to volume, it is also a critical conversion factor for estimating SOC stocks. Even small methodological biases in BD determination can therefore propagate into substantial uncertainties in SOC stock estimates, with direct implications for soil health assessments, carbon budgeting, and evaluations of land management and climate mitigation strategies. Despite its critical importance, BD is frequently one of the most commonly missing variables in standard soil datasets. BD can be determined using a range of direct and indirect methods, each involving trade-offs among accuracy, cost, and operational feasibility. Direct approaches (e.g. core, clod, or excavation methods) are widely applied but are labour-intensive and sensitive to operator technique, sampling depth, cylinder dimensions, and soil moisture conditions. Indirect approaches, including pedotransfer functions (PTFs), can reduce field effort and, in some cases, improve spatial coverage, but they require careful calibration and high-quality ancillary data (e.g. texture, organic carbon, climatic variables) and may introduce depth-dependent uncertainties. These issues are particularly critical in heterogeneous Mediterranean landscapes, where BD exhibits strong spatial and vertical variability.
Within this context, a legacy data-rescue activity was conducted for Sicily, the largest island of the Mediterranean Basin, integrating multiple legacy datasets and the first results from the Soils4MED survey, to compile approximately 2,000 records of soil texture, SOC, and BD. Data were collected from the soil province geodatabase of Italy, the Land Use and Coverage Area frame Survey (LUCAS) topsoil, the European Hydropedological Data Inventory (EU-HYDI), the Soils4MED project, and several peer-reviewed papers. To ensure consistency across heterogeneous sampling depths, data harmonization was performed using a mass-preserving approach based on the equal-area quadratic spline method.
To generate spatially explicit estimates, a geostatistical approach was used to create maps at 1 km and 500 m resolution of all input variables of the PTFs before the determination of BD to reduce the propagation of errors in input data. These products, together with the original point data, were compared against established European-scale BD and packing density datasets to identify potentially biased areas. Moreover, an uncertainty BD map was generated by a geostatistical stochastic simulation to provide the quality of BD assessments at the different locations of Sicily. The results highlight the added value of data integration, geostatistics, and PTFs for improving BD representation and supporting robust SOC stock assessments in Mediterranean soils, in line with the objectives of the EU Soil Mission.
Acknowledgements
The SOILS4MED project is part of the PRIMA programme supported by the European Union.