Modeling Soil Organic Carbon Changes Using Signal-To-Noise Analysis: A Case Study Using European Soil Survey Datasets.

Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Glob Chang Biol . 2026 Apr 10;32(4):e70813. doi: 10.1111/gcb.70813 Search in PMC Search in PubMed View in NLM Catalog Add to search Modeling Soil Organic Carbon Changes Using Signal‐To‐Noise Analysis: A Case Study Using European Soil Survey Datasets Xuemeng Tian Xuemeng Tian 1 OpenGeoHub, Doorwerth, the Netherlands 2 Laboratory of Geo‐Information Science and Remote Sensing, Wageningen University & Research, Wageningen, the Netherlands Find articles by Xuemeng Tian 1, 2, ✉ , Sytze de Bruin Sytze de Bruin 2 Laboratory of Geo‐Information Science and Remote Sensing, Wageningen University & Research, Wageningen, the Netherlands Find articles by Sytze de Bruin 2 , Florian Schneider Florian Schneider 3 Thünen Institute of Climate‐Smart Agriculture, Germany Find articles by Florian Schneider 3 , Martin Herold Martin Herold 2 Laboratory of Geo‐Information Science and Remote Sensing, Wageningen University & Research, Wageningen, the Netherlands 4 Helmholtz GFZ German Research Centre for Geosciences, Remote Sensing and Geoinformatics, Potsdam, Germany Find articles by Martin Herold 2, 4 , Kirsten de Beurs Kirsten de Beurs 2 Laboratory of Geo‐Information Science and Remote Sensing, Wageningen University & Research, Wageningen, the Netherlands Find articles by Kirsten de Beurs 2 Author information Article notes Copyright and License information 1 OpenGeoHub, Doorwerth, the Netherlands 2 Laboratory of Geo‐Information Science and Remote Sensing, Wageningen University & Research, Wageningen, the Netherlands 3 Thünen Institute of Climate‐Smart Agriculture, Germany 4 Helmholtz GFZ German Research Centre for Geosciences, Remote Sensing and Geoinformatics, Potsdam, Germany * Correspondence: Xuemeng Tian ( [email protected] ) ✉ Corresponding author. Revised 2026 Jan 23; Received 2025 Aug 4; Accepted 2026 Jan 26; Issue date 2026 Apr. © 2026 The Author(s). Global Change Biology published by John Wiley & Sons Ltd. This is an open access article under the terms of the http://creativecommons.org/licenses/by/4.0/ License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. PMC Copyright notice PMCID: PMC13066779  PMID: 41958231 ABSTRACT Soil organic carbon (SOC) is a key indicator of soil health and a crucial component of climate mitigation, making its reliable monitoring increasingly important. While Digital Soil Mapping (DSM) based on Machine Learning and Earth Observation (EO) data enables the generation of time series of spatially explicit SOC predictions, detecting temporal changes from these model predictions remains challenging due to the relatively large associated uncertainties. Although prediction uncertainties are now commonly reported, few studies have explicitly accounted for them when assessing SOC change. This study introduces a model‐based signal‐to‐noise ratio (SNR) framework to assess the detectability of SOC change using both the state‐first approach—modeling SOC states at each time point and then deriving change—and the change‐first approach—modeling SOC change directly from repeated measurements. SNR is defined as the ratio of predicted SOC (concentration, g/kg) change to its modeled uncertainty, enabling evaluation of change‐model reliability at pixel levels. Applied to repeated SOC observations from the pan‐European Land Use and Coverage Area Frame Survey, this framework assesses the reliability of SOC change modeling across multiple land‐cover types using Random Forest and Quantile Regression Forests. At the site level, prediction accuracy was poor and SNR values were consistently low. An illustrative aggregation analysis showed that spatial averaging improved SNR, supporting SOC change assessments at broader scales. However, further work is needed to incorporate land use and management information and to systematically examine how different aggregation schemes affect the results in various contexts, ensuring that aggregated outcomes remain meaningful and policy‐relevant. As an internal metric based on model predictions and their estimated uncertainty, SNR provides a practical diagnostic of change‐model confidence, especially when repeated ground‐truth SOC measurements are not available. We advocate for routine SNR reporting to enhance the transparency and credibility of DSM‐based SOC change monitoring. Keywords: change detection, digital soil mapping, machine learning, signal‐to‐noise ratio, soil organic carbon This conceptual illustration underpins the Signal‐to‐Noise Ratio (SNR) framework proposed in the study to assess SOC change detectability using repeated Soil Organic Carbon (SOC) observations, Machine Learning and Earth Observation data. Using a simple simulated time series, the figure summarizes the two modelling approaches evaluated in this study: State‐first approach: Models are trained to predict SOC states at each time step. SOC change is then derived by differencing the predicted states, and state‐level uncertainties are propagated to obtain the uncertainty of the derived change. Change‐first approach: SOC change between measurements is first computed, and ML models are trained to predict the change directly along with its associated uncertainty. By treating the magnitude of the estimated change as the signal and the corresponding uncertainty as the noise, SNR can be computed for both approaches. SNR indicates how confidently SOC change can be detected: high SNR reflects a reliable change signal, whereas low SNR indicates that uncertainty overwhelms the predicted change. 1. Introduction Soil organic carbon (SOC) is widely recognized as a key indicator of soil health and a critical component of climate mitigation, owing to its role in sequestering atmospheric carbon (Lal 2004 ; Lehmann et al. 2020 ; Stockmann et al. 2013 ; Don et al. 2024 ). Reliable information on SOC dynamics is therefore increasingly sought by scientists, policymakers, and land managers, particularly for monitoring, reporting, and verification (MRV) purposes (Lehmann et al. 2020 ; Smith et al. 2020 ). To meet this growing need for reliable SOC information, digital soil mapping (DSM) has gained momentum in recent years, particularly through approaches combining machine learning (ML) and satellite‐based Earth observation (EO) data (Wadoux et al. 2020 ). ML captures complex, nonlinear relationships among predictors, while EO provides consistent, high‐resolution, and regularly updated surface information. Together with the expansion of soil survey databases, these advances have driven a surge in spatiotemporal SOC mapping (Hengl et al. 2018 ; Wadoux et al. 2020 ). Notable examples include the biannual Worldsoils global maps (van Wesemael et al. 2024 ), the SoilHealthDataCube for Europe (Tian, de Bruin, et al. 2025 ), and national‐scale time series for Argentina, the Netherlands, South Africa, Hungary, and others (Heuvelink et al. 2021 ; Helfenstein et al. 2024 ; Venter et al. 2021 ; Szatmári, Laborczi, et al. 2024 ). Beyond reconstructing past SOC dynamics, several studies have also extended similar approaches to project future trends (Yigini and Panagos 2016 ; Wang et al. 2022 ; Zhang et al. 2023 ; Soils Revealed 2025 ). From such SOC maps, SOC change is derived by first predicting SOC levels at each time step and then calculating their difference or fitting a trend line (Chen et al. 2019 ; Yang et al. 2022 ; Li et al. 2022 ; Helfenstein et al. 2024 ; Szatmári, Pásztor, et al. 2024 ; Meng et al. 2024 ), an approach here referred to as the state‐first strategy. In contrast, the change‐first approach derives SOC changes directly from repeated observations and models those changes as the target variable. By focusing on change rather than absolute states, this strategy can mitigate correlated noise between time steps (De Rosa et al. 2024 ). It remains less common due to the scarcity of long‐term, repeated SOC measurements. Practical mapping applications are still rare, with De Rosa et al. ( 2024 )'s analysis of repeated Land Use and Coverage Area Frame Survey (LUCAS, Orgiazzi et al. 2018 ) topsoil samples being one of the few examples. Uncertainty in SOC prediction can arise from multiple sources, including soil sampling variability, EO data quality, and model assumptions (Heuvelink et al. 2021 ; Zhang, Huang, and Yang 2024 ; Even et al. 2025 ). As a result, accounting for prediction uncertainty has become an integral part of DSM (G. B. Heuvelink 2014 , 2018 ). Various methods have been developed to quantify model uncertainty, among which quantile regression forests (QRF) are widely applied (Vaysse and Lagacherie 2017 ; Wadoux et al. 2020 ). Consequently, most recent DSM products now report uncertainty alongside predicted values, typically as prediction intervals derived from model distributions (G. B. Heuvelink 2014 ; Kasraei et al. 2021 ; Schmidinger and Heuvelink 2023 ). Although uncertainty is particularly relevant when expected SOC changes are small relative to prediction uncertainty (Poeplau and Don 2013 ; Paustian et al. 2019 ), few studies have quantitatively tested whether modeled SOC changes exceed their associated uncertainty. In the state‐first approach, uncertainty is often visualized through prediction intervals alongside SOC time series (Tian, de Bruin, et al. 2025 ), which provides useful context but remains qualitative. A rare quantitative example is Szatmári, Pásztor, et al. ( 2024 ), who considered uncertainty in SOC change but focused only on trend direction (positive vs. negative). For the change‐first approach, applications are even rarer; in De Rosa et al. ( 2024 )'s study, uncertainty was shown alongside predicted SOC change, but only at an aggregated longitudinal level. To address this gap, we introduce a model‐based signal‐to‐noise ratio (SNR) framework to assess whether predicted SOC changes are distinguishable from model uncertainty. SNR is defined as the ratio between the magnitude of the predicted SOC change (the signal) and its modeled prediction uncertainty (the noise). A high SNR (SNR > 1) indicates that the modeled change exceeds its uncertainty, providing a straightforward internal measure of model confidence—particularly valuable when independent repeated soil samples for external validation are scarce. We demonstrate the approach with repeated SOC measurements from the pan‐European LUCAS topsoil (0–20 cm) dataset, which enables evaluation across diverse land‐cover types under both the state‐first and change‐first approaches. Although implemented at the pixel level, the SNR framework can also be applied to aggregated spatial units (e.g., administrative boundaries or regular grids), which is particularly relevant for MRV applications at different spatial levels. Spatial aggregation further improves SNR by dampening localized noise (Wadoux and Heuvelink 2023 ; Szatmári, Pásztor, et al. 2024 ; Tian, de Bruin, et al. 2025 ), a point we demonstrate through case studies in this paper. Although SOC stock (e.g., Mg/ha) is more relevant for climate change mitigation, repeated measurements of bulk density, coarse fragments and depth to bedrock are even scarcer than those of SOC concentration (e.g., g/kg). Hence, this study focuses on SOC concentration, with “SOC” hereafter referring to concentration. While SOC concentration is used in this study due to current data limitations, the proposed SNR framework can, in principle, be extended to SOC stock when the necessary auxiliary data become available. Through this SNR‐based evaluation, we address the following research questions, which together aim to assess the reliability and detectability of SOC change modeling across spatial scales: RQ1 : How do SOC change prediction capabilities across Europe compare to those reported for the Netherlands (Helfenstein et al. 2024 ) and German croplands (Broeg et al. 2024 )? RQ2 : How do the change‐first and state‐first strategies differ in predictive accuracy and SNR when modeling SOC change? RQ3 : Which factors influence SNR in SOC change modeling, such as land cover, monitoring duration, sampling frequency, and spatial scales of interest? RQ4 : How does SNR relate to prediction accuracy in SOC change modeling? 2. Materials and Methods 2.1. Definitions and Notation for SOC and Its Change We denote the measured SOC (g/kg) at location s and time t i as c s , t i . SOC change can be described in two main ways: (i) as the net difference between paired observations over a defined interval ( δ s , t 1 → t 2 ; g/kg) (Lettens et al. 2005 ; Poeplau and Don 2013 ; Harbo et al. 2023 ), and (ii) as a temporal trend estimated through regression models that relate SOC to time. In the absence of detailed knowledge about temporal SOC dynamics, it is commonly assumed that SOC changes linearly with time, and linear models are fitted across multiple time points. The resulting slope parameters, denoted β s , t 1 → t n (g/kg/year) (Bellamy et al. 2005 ; Harbo et al. 2025 ) represent the average rate of SOC change over the time window t 1 → t n for spatial unit s . All SOC change metrics are derived from repeated measurements at the same location ( c s , t i ). For simplicity, spatial and temporal indices (e.g., s , t ) are omitted throughout the manuscript unless explicitly required. Accordingly, SOC measurements and their derived metrics are denoted as c , δ , and β , representing the observed SOC concentration, net change, and temporal trend, respectively. Their corresponding model‐predicted values are denoted as c ^ , δ ^ , and β ^ . 2.2. Data Preparation In this study, we primarily used soil data from the LUCAS soil database, supplemented with three national soil datasets. Two of these are from Spain: Parcelas COS (top 30 cm depth) and Parcelas INES (top 10 cm) (Serrano et al. 2022 ). The third dataset, BZE‐LW, is from Germany and represents the core dataset of the first German agricultural soil inventory (Poeplau et al. 2020 ), containing measurements across multiple soil depths. As the only database providing repeated topsoil samples (0–20 cm) with pan‐European coverage, LUCAS served as the reference dataset and foundation of this study. The analysis focused on topsoil samples within the 0–20 cm depth range. For sites with measurements at multiple depths, SOC values were interpolated to a standardized depth of 10 cm—representing the midpoint of the target layer. To preserve the shape of data and avoid overshooting, we used the Piecewise Cubic Hermite Interpolating Polynomial method (PCHIP) (Fritsch and Butland 1984 ), which achieved a high coefficient of determination (0.8) in reconstructing SOC values at missing depths in an internal accuracy assessment ( Supplementary Notebook 02e ). The final harmonized dataset includes 87,684 c observations, from all over Europe (Figure 1A ), covering year 2002 to 2019 (Figure 1C ). Summary statistics are provided in Table 1 . FIGURE 1. Open in a new tab Spatial density of all point measurements (A) and filtered time‐series observations (B) aggregated to 50 km grid cells across Europe, together with their temporal distribution across survey years (C). TABLE 1. Summary statistics of point SOC (g/kg) per dataset. Source Year N sites N samples Mean Median SD LUCAS 2009, 2012, 2015, 2018 27,737 62,437 45.98 20.70 82.03 Parcelas INES 2002–2018 21,502 21,502 36.63 23.53 38.67 BZE‐LW 2011–2018 2953 2953 30.07 18.37 49.22 Parcelas COS 2019–2020 789 789 40.78 27.17 34.64 Open in a new tab Among the LUCAS observations, a subset includes repeated measurements taken at the same locations over three sampling rounds: 2009/2012, 2015, and 2018, enabling the analysis of SOC dynamics over time. To ensure data consistency and quality, a filtering procedure was applied to keep only consistent SOC time series. The following criteria were used: (1) each time series must contain at least three observations, representing the longest possible duration; (2) the maximum spatial distance between repeated measurements must not exceed 30 m, a practical threshold aligning with the resolution of our previous SOC maps (Tian, de Bruin, et al. 2025 ), under the simplifying assumption that points within 30 m fall within the same map pixel; and (3) the standard deviation of the carbon‐to‐nitrogen (C/N) ratio across the time series must be less than 3.8, based on the observed bimodal distribution of C/N differences, where 3.8 marks a local minimum between two peaks (see Supplementary Notebook 02b for details; Ziche et al. 2022 ). This filtering yielded 10,204 SOC time series from the LUCAS dataset, covering most of Europe (Figure 1B ), with observations primarily from the years 2009, 2015, and 2018 (Figure 1C ). A small subset of time series includes measurements from 2009, 2012, and 2018, though this group is not clearly visible due to its limited size. The change ( δ ) was calculated for all valid pairs within each time series, including step‐to‐step (e.g., t 1 → t 2 , t 2 → t 3 ) and longer‐interval (e.g., t 1 → t 3 ) differences. This resulted in a total of 30,612 δ values across the dataset. We estimated the rate of change at each location using the Theil–Sen estimator—a non‐parametric, robust linear regression method well‐suited for small sample sizes and resistant to outliers (Theil 1950 ; Sen 1968 ). Slopes were derived from the full SOC time series at each site, yielding 10,204 β values. This was implemented using TheilSenRegressor from sklearn.linear_model . An independent test set of 1800 SOC time series (about 20% of the dataset, corresponding to 5400 point observations and 5400 time‐step pairs) was randomly selected from the filtered LUCAS dataset. This size was chosen to ensure a sufficiently large and diverse test set for evaluation, while leaving enough data for training and exploration. This test set was used for evaluating prediction accuracy, assessing uncertainty estimates, and conducting SNR analysis. 2.3. Covariates A wide range of covariates (also referred to as predictors or features ) was included to represent environmental factors influencing soil formation, as described by the SCORPAN model. Specifically, SCORPAN refers to s: soil properties; c: climate; o: organisms (vegetation); r: topography; p: parent material (lithology); a: age; and n: spatial position (McBratney et al. 2003 ). Topographic variables (Ho et al. 2025 ) and parent lithology (Isik et al. 2024 ) were included to represent terrain‐driven controls on hydrology, erosion, and soil formation processes. Plant functional type (Harper et al. 2023 ) and vegetation cover fraction (Sun et al. 2023 ) maps were used to characterize vegetation composition and cover, which are directly related to carbon inputs to the soil (Smith 2008 ). Climate variables include precipitation (Karger et al. 2021 ), temperature (Wan 2006 ), and water vapor (Lyapustin and Wang 2018 ), as well as bioclimatic variables (Karger et al. 2017 ) that represent biologically meaningful climate constraints relevant to vegetation growth and soil carbon processes. Thematic peatland (Widyastuti et al. 2024 ) and cropland (Potapov et al. 2022 ) datasets were included to represent persistently waterlogged peat soils with high SOC contents and intensively managed croplands with direct human intervention. Soil moisture (Bauer‐Marschallinger et al. 2018 ) was included to indicate water availability and aeration conditions that regulate organic matter decomposition and SOC stabilization. In addition, less‐processed remote sensing signals were incorporated, including Landsat reflectance and derived spectral indices (Tian, Consoli, et al. 2025 ) and bare‐soil composites from Sentinel imagery (Rogge et al. 2018 ). Two Synthetic Aperture Radar (SAR) backscatter signals at different wavelengths (Shimada and Ohtaki 2010 ; Wagner et al. 2021 ) were also included. These signals are sensitive to surface properties and vegetation structure and have been shown to be informative for modeling SOC in superficial soil layers (Baghdadi et al. 2008 ). Some covariates were available as time series, providing multiple observations for the same location over time, including vegetation‐related variables, climate variables, cropland information, and soil moisture. Other covariates provided only a single value per location, including static variables or long‐term summaries. Static layers include slowly varying properties such as topography and lithology, as well as variables treated as static as a pragmatic compromise given data availability and consistency at the continental scale, including Synthetic Aperture Radar (SAR) backscatter and peatland cover. Long‐term summaries were derived from multi‐year time series for selected covariates and include statistical descriptors such as the mean, percentiles, and standard deviation, resulting in a single value per location. In this study, “long‐term” refers to the period from 2000 to 2022, corresponding to the time span over which most dynamic covariates are consistently available. The resulting initial feature pool consisted of 582 individual covariates, organized into 15 covariate groups, which in turn represent seven broader SCORPAN classes (Table 2 ). The full list of covariates is not included in the main text but is available in the Supporting Information at our Supplemetary Repository . Various preprocessing steps were applied to the original datasets, including temporal aggregation, gap‐filling of missing values, and reprojection to ensure consistent coordinate reference systems, spatial coverage, and temporal alignment across all layers. TABLE 2. Covariate groups used in this study. Class Covariate group Source Temporal resolution Spatial resolution Topography Digital terrain model and derived land surface parameters Ho et al. ( 2025 ) Static , 60, 120, 240, 480, and 960 m Lithology Lithology type probability Isik et al. ( 2024 ) Static m Vegetation Plant functional type Harper et al. ( 2023 ) Annual m Vegetation Vegetation cover fraction Sun et al. ( 2023 ) Annual m Climate Precipitation Karger et al. ( 2021 ) Annual km Climate Land surface temperature Wan ( 2006 ) Annual km Climate Water vapor Lyapustin and Wang ( 2018 ) Annual km Climate Bioclimate Karger et al. ( 2017 ) Long‐term km Imagery Landsat spectral index Tian, Consoli, et al. ( 2025 ) Annual, long‐term m Imagery Bare‐soil composite Rogge et al. ( 2018 ) Long‐term m Imagery Sentinel‐1 SAR backscatter Wagner et al. ( 2021 ) Static m Imagery PALSAR backscatter Shimada and Ohtaki ( 2010 ) Static m Land cover Peatland Widyastuti et al. ( 2024 ) Static km Land cover Cropland Potapov et al. ( 2022 ) Annual m Soil Soil moisture Bauer‐Marschallinger et al. ( 2018 ) Annual km Open in a new tab Note: Note that each entry represents a group of covariate layers, which may include multiple datasets with different temporal and spatial resolutions. These covariate layers were overlaid with c observations based on spatial coordinates (latitude and longitude) and temporal information (year) to construct the training matrix for the R F c . To reduce dimensionality and improve interpretability, we applied the Repeated Subsampling‐Based Cumulative Feature Importance (RSCFI) method (Tian, de Bruin, et al. 2025 ). This approach, conceptually similar to recursive feature elimination with cross‐validation, iteratively removes features with low cumulative importance scores, thereby retaining only the most relevant predictors while balancing model performance and computational efficiency. Details of the RSCFI implementation are provided in Supplementary Notebook 07a . The covariates identified as important by R F c were then organized to reflect environmental changes over time, supporting the construction of training matrix for R F δ and R F β . Following Yang et al. ( 2022 ), we applied time‐series feature engineering to convert state‐based covariates into change‐based representations. Each covariate time series was partitioned into three periods: (i) the change interval (between start and end years), (ii) the years preceding the change (from 2000—the earliest common year—up to the start year), and (iii) the full period from 2000 to the end year. From each period, we derived summary statistics of feature values—mean, standard deviation, and linear slope (as a proxy for trend)—yielding a structured set of dynamic predictors aligned with the temporal spans of δ and β . These were then used to train R F δ and R F β . 2.4. SOC Model Approaches All models are implemented using RF and QRF, two widely used methods in DSM for SOC prediction and uncertainty quantification (Hengl et al. 2018 ; Wadoux et al. 2020 ). We used a modified version of the RF algorithm to capture full prediction distributions by storing values from all valid nodes across trees. This mimics the behavior of QRF, enabling simultaneous estimation of the target variable (via the mean) and associated uncertainty (via the standard deviation). This was implemented by modifying the RandomForestRegressor in scikit‐learn to keep node‐level outputs. Three separate RF models were constructed using the same workflow, including feature selection, hyperparameter tuning, and model fitting, to predict c , δ and β , respectively. state‐first : This approach relies on a single model, R F c , to generate spatial–temporal c ^ , from which δ ^ and β ^ are subsequently derived. – Model of spatial–temporal point SOC values ( R F c ) : The training matrix was constructed by aligning SOC observations with the corresponding environmental features, drawn from 15 covariate groups (Table 2 ), based on spatial coordinates (latitude and longitude) and temporal information (year). change‐first : In this approach, the target variable directly corresponds to the modeled output. Temporal dynamics are incorporated through feature engineering of covariates, aiming to reflect underlying processes rather than static states (see next section for details). – Model of paired SOC change ( R F δ ) : The training matrix was constructed by matching δ values with features that reflect environmental changes over the corresponding time period at the same location. – Model of SOC trend slope ( R F β ) : As with the R F δ , covariates were engineered to reflect environmental processes across the time span of each series. Performance metrics—including mean absolute error (MAE), bias, and concordance correlation coefficient (CCC)—were calculated on an independent test set to evaluate the prediction accuracy of the three models: R F c , R F δ , and R F β , used under both the state‐first and change‐first approaches. This resulted in five prediction scenarios (Table 3 ). TABLE 3. Overview of RF models used in the state‐first and change‐first approaches and their corresponding target variables. Approach Model Target Predicted State‐first R F c c c ^ δ δ ^ derived as temporal differences of c ^ β β ^ derived as temporal slopes fitted to c ^ Change‐first R F δ δ δ ^ R F β β β ^ Open in a new tab We used the standard deviation of RF node outputs as a proxy for prediction noise in the SNR calculation—an approach valid only if this spread reliably reflects uncertainty. Therefore, the models were also evaluated for the validity of their uncertainty estimates using two metrics: the Prediction Interval Coverage Probability (PICP), which measures the proportion of observed values falling within the predicted intervals; and the Quantile Coverage Probability (QCP), which assesses the symmetry of the uncertainty distribution by evaluating single‐quantile predictions. Accuracy plots of PICP and QCP—scatter plots comparing observed and expected coverage rates across a range of expected prediction intervals—were used to visualize and evaluate local uncertainty estimation performance, as recommended by Schmidinger and Heuvelink ( 2023 ). 2.5. SNR Analysis The SNR was defined as the ratio between the model's predicted value and its associated prediction uncertainty. SNR thus quantifies the extent to which a modeled change exceeds its own uncertainty, providing an internal diagnostic of model confidence. Its specific computation varies by model, target variable, and approach, as shown in Table 4 , where the numerator represents the signal —the absolute magnitude of the predicted target ( ∣ c ^ ∣ , ∣ δ ^ ∣ , or ∣ β ^ ∣ ). Although the calculation of SNR varies across models, all rely on the prediction distribution—that is, the set of terminal node values from the QRF. This conditional distribution is assumed to approximate model uncertainty within the local feature space. TABLE 4. SNR calculation formulas for each model and target variable. Target Model R F c R F δ R F β R = c ^ 1 , c ^ 2 , … , c ^ n R = δ ^ 1 δ ^ 2 … δ ^ n R = β ^ 1 β ^ 2 … β ^ n c c ^ std R — — δ ∣ δ c ^ t 0 c ^ t 1 ∣ std R t 0 2 + std R t 1 2 * ∣ δ ^ ∣ std R # — β ∣ β c ^ t 0 c ^ t 1 c ^ t 2 ∣ std β R t 0 R t 1 R t 2 * — ∣ β ^ ∣ std R # Open in a new tab Note: R denotes the prediction distribution, that is, all prediction instances for a target variable obtained from the full valid terminal node values of the RF. * Indicates the state‐first approach. # Indicates the change‐first approach. In the change‐first approach, noise is quantified as the standard deviation of this prediction distribution. While this is similar to QRF, which uses the same prediction set to estimate conditional quantiles (Meinshausen 2006 ), our method focuses on capturing the uncertainty level via spread. This procedure is applied consistently across the three direct models ( R F c , R F δ , and R F β ), corresponding to the diagonal entries in Table 4 . For the state‐first approach, the R F c first predicts c ^ values, from which both δ ^ and β ^ are derived (Table 3 ). To estimate noise in δ ^ , we assumed the two c ^ estimates to be independent, as RF models treat samples independently. Under this assumption, the noise for δ ^ is computed as the square root of the sum of variances from both time steps. To estimate noise for β ^ , we implemented a Monte Carlo simulation. For each iteration, c ^ values were randomly sampled from the full prediction distribution (i.e., all valid leaf node outputs) at each time step, and a linear regression was fitted to estimate the corresponding slope. This process was repeated 200 times, generating a distribution of β ^ values. The standard deviation of this distribution was used as the noise estimate for β ^ . To explore the relationship between internal model confidence (SNR) and actual predictive accuracy, we compared SNR with relative absolute error (RAE) for each δ ^ and β ^ in the validation dataset. RAE was used to evaluate prediction accuracy across all target variables ( c , δ , and β ). It is defined as the absolute difference between the predicted and observed values, normalized by the observed value. 2.6. Spatial Aggregation While spatial averages can be calculated by averaging predictions over the area of interest (AOI), quantifying the corresponding uncertainty is less straightforward—particularly owing to the spatial dependence of prediction errors. In this study, we adopted an approach similar to that used by Araza et al. ( 2022 ) and Wadoux and Heuvelink ( 2023 ), which estimates the uncertainty of spatial aggregates using model‐based point‐support uncertainty‐prediction standard deviation, in our case, derived from discretization points within the AOI. This method avoids the numerical complexity of traditional statistical techniques, making it especially suitable for large‐scale applications. Specifically, it is implemented as follows: sd AOI = 1 B 2 ∫ s ∈ B ∫ u ∈ B sd s · sd u · ρ s − u ds du (1) The B in Equation ( 1 ) represents the number of discretization points, generated from the AOI; sd s and sd u denote the standard deviations of the point support prediction residuals at locations s and u , assuming the prediction errors are proportional to the prediction uncertainty (standard deviation of prediction distribution from QRF) at point support. The ρ s − u represents the correlation function of the standardized prediction error at the separation distance between s and u . Standardized prediction errors and corresponding uncertainty were derived through cross‐validation, where predictions were made on samples excluded from model training. Errors were calculated as the difference between observed and predicted values, then standardized by dividing by the prediction uncertainty. We assumed second‐order stationarity of the standardized errors—that is, their spatial correlation depends only on the Euclidean distance between locations. The correlation structure was estimated by fitting a variogram γ h using: ρ h = sill − γ h sill , (2) where sill denotes the fitted sill of the variogram. To assess the effect of spatial aggregation on SNR, SOC prediction maps at 1 km 2 resolution were first generated, which served as the basis for subsequent aggregation. We considered four levels of spatial support, a concept from geostatistics describing the area over which a variable is aggregated: 2.5, 10, 40, and 200 km 2 . These scales, spanning fine to coarse resolutions, were chosen to demonstrate the effect of spatial aggregation on SNR and are not intended as prescriptive values. The analysis focused on Germany and Spain, where additional data is available. Each country was divided into square grid units corresponding to the four support sizes. Within each grid unit, map pixels were used as discretization points for aggregation, and B denotes the number of pixels contained in each grid unit at the respective spatial support. Spatial aggregation was applied only to the R F c model under the state‐first approach. This restriction was mainly due to the limited density of available SOC time series measurements, which must be sufficient to fit standardized error variograms for aggregated uncertainty estimation. Moreover, as state‐first approach remains the standard in SOC change modeling, results are broadly applicable to current SOC time series products. 3. Results 3.1. Prediction Accuracy and Uncertainty Estimation Evaluation The prediction accuracy of c ^ was first evaluated to justify the use of the state‐first approach. R F c achieves a CCC of 0.344 in the original scale (Figure 2A ) and 0.705 in the log‐transformed scale (Figure 2B ). These values are comparable to those reported in previously published continental‐scale DSM models (Tian, de Bruin, et al. 2025 ; van Wesemael et al. 2024 ). While the model shows a slight positive bias—tending to overestimate SOC—it consistently underestimates the highest observed SOC values. FIGURE 2. Open in a new tab Prediction accuracy of c using R F c . Panel (A) shows the density plot of observed versus predicted values on the original scale (g/kg), and panel (B) shows the same results with log‐transformed values (log1p transformed g/kg) for greater detail at lower concentrations. Both panels are based on the same validation dataset. Performance metrics are calculated and presented in each panel. Figure 3 shows the accuracy of SOC change predictions— δ ^ (top) and β ^ (bottom)—using both state‐first (left) and change‐first (right) approaches. Unlike the good performance for c ^ , δ ^ and β ^ show poor accuracy across all models and approaches, with all models yielding near‐zero CCC and high MAE relative to the scale of the target variables. For state‐first (Figure 3A,C ), predictions are weakly correlated with observations. Accuracy does not improve under change‐first (Figure 3B,D ). Both approaches struggle to reliably capture the correct direction of change. The state‐first and change‐first approaches yield conflicting biases for both δ ^ and β ^ . FIGURE 3. Open in a new tab Prediction accuracy of δ and β using two approaches: State‐first based on predictions from R F c (A and C), and change‐first using R F δ (B) and R F β (D). Performance metrics are shown in the upper‐left corner of each subplot. Despite overall poor performance, state‐first approach (Figure 3A ) shows a wider spread of δ ^ predictions than change‐first (Figure 3B ), where predictions cluster around zero. This contrast is less pronounced for β ^ (Figure 3C,D ). All three models demonstrated reasonable performance in estimating uncertainty when applied to their direct target under the change‐first approach (Figure 4 ). A slight underestimation—indicative of over‐optimism—was observed in the PICP of the R F δ and R F β (Figure 4D,F ) when the expected prediction interval was around 50%. When assessing the sharpness of prediction intervals, the results (see Supplementary Notebook 15 ) show that the uncertainties estimated by R F δ and R F β under the change‐first approach are generally narrower but exhibit less variation. FIGURE 4. Open in a new tab Accuracy plots of QCP (left column, A, C and E) and PICP (right column, B, D and F) used to evaluate the R F c (A and B), R F δ (C and D), and R F β (E and F)—in estimating uncertainty for their respective target variables. All plots present scatter points comparing observed coverage rates on the y‐axis with expected coverage rates on the x‐axis across a range of prediction intervals. 3.2. SNR Analysis of SOC and Its Change Modeling Model‐based SNR provides an internal diagnostic of SOC change detectability, especially when repeated ground‐truth data are unavailable. The SNR analysis for c ^ is presented in Figure 5 , stratified by land‐cover type. All Level 1 land cover classes from the LUCAS survey are represented in the study, except for Wetland and Water Areas and Artificial Land , which were excluded due to insufficient sample sizes (fewer than 10 observations in the test set). As shown in Figure 5 , the histograms of observed and predicted SOC generally align for Bareland , Cropland , Grassland , and Shrubland , while a clear mismatch is observed for Woodland . In Woodland , the predicted SOC values are shifted toward higher concentrations and appear compressed, with less extreme values at both ends. FIGURE 5. Open in a new tab Analysis of the SNR for c ^ by R F c . Left: Histograms of observed and predicted SOC values across different land‐cover types. Right: Summary statistics of sample count, signal ( c ^ ), error, noise, and SNR derived from the independent test set for each land‐cover type. In the right panel, boxes indicate the interquartile range, and whiskers represent the remaining distribution excluding outliers. Note that both plots are based solely on the independent test set. As the independent test set is a random subset of the full dataset, the sample counts reflect the original data distribution, with the majority of observations coming from Cropland , followed by Woodland and Grassland . Bareland and Shrubland are minimally represented. The signal component increases from Bareland to Woodland , consistent with increasing vegetation content. Both prediction errors and noise estimates also increase along this gradient. In terms of SNR, Cropland shows the highest values, while Woodland has the lowest, despite having the second‐highest number of available samples. SNR values are mostly above one for Bareland and Cropland , close to one for Grassland and Shrubland , and predominantly below one for Woodland . Figure 6 presents the SNR of δ ^ (bottom panel) for both the state‐first and change‐first approaches, along with the corresponding signal, noise, and the distribution of observed δ for reference. Observed δ values tend to be larger and more variable in land‐cover types with higher vegetation. For example, δ is larger when Woodland is involved—either when the land cover remained as Woodland across both sampling years or transitioned to or from it. In contrast, transitions involving Cropland show lower magnitudes and narrower spreads. Although the exact values of δ and δ ^ do not match at the individual prediction level (Figure 3 ), both the state‐first and change‐first approaches reproduce the overall trend of δ in δ ^ across land‐cover types, with higher predicted signal magnitudes for transitions involving vegetation‐rich classes. FIGURE 6. Open in a new tab Summary statistics of observed δ , as well as signal, noise, and SNR of δ ^ , comparing results from the R F δ under the change‐first approach and the R F c under the state‐first approach, grouped by land cover transition pairs. Bars represent median values, and whiskers indicate the 5th to 95th percentile range. Land cover codes are as follows: B = Bareland, C = Cropland, G = Grassland, S = Shrubland, and W = Woodland. For example, “CC” indicates that the location was classified as Cropland in both sampling years. Based on the general increase in vegetation content along the gradient B → C → G → S → W, land cover transitions are grouped into three categories shown in the top panel: No change, increasing vegetation content, and decreasing vegetation content. Overall, the change‐first approach yields lower SNR values than the state‐first approach, except for transitions involving Woodland (e.g., WW, SW, and GW). The variation in SNR across land‐cover types is more pronounced for the change‐first approach. This approach reproduces the signal pattern across land‐cover classes—where vegetation‐rich covers correspond to larger signals—while maintaining relatively stable noise estimates. Consequently, transitions involving high‐vegetation land covers show higher SNR values and greater variation across classes. In contrast, under the state‐first approach, although the strongest and most variable signals also occur in high‐vegetation classes such as Woodland , the noise magnitudes follow similar land‐cover patterns as the signal and δ . As a result, SNR values remain relatively consistent across transitions and generally below 0.4. Compared to δ ^ , the SNR estimates for β ^ from both approaches show stronger agreement, with more similar magnitudes and variation across land‐cover types, although the overall SNR patterns for β remain similar to those for δ (Figure 7 ). Observed β values are higher and more variable for transition types involving vegetation‐rich land covers, a pattern mirrored in the predicted signals from both the state‐first and change‐first approaches. For the change‐first approach, noise estimates remain relatively stable, resulting in larger SNR variations across land‐cover transitions compared with the state‐first approach. FIGURE 7. Open in a new tab Summary statistics of signal (top panel), noise (middle panel), and SNR (bottom panel) for β , comparing results from the β ‐model under the change‐first approach and the c ‐model under the state‐first approach, grouped by land cover time series. Bars represent median values, and whiskers indicate the 5th to 95th percentile range. Land cover codes follow those used in Figure 6 , and time series are grouped into three categories: No change, increasing vegetation content, and decreasing vegetation content. For both δ ^ and β ^ , only a weak positive correlation was observed between SNR and RAE (Spearman's rank correlation coefficient was 0.17 in both cases for the state‐first approach and close to zero for the change‐first approach; see Supplementary Notebook 12 ). The majority of predictions fall into the “high RAE and low SNR” quadrant, indicating that large prediction errors are generally associated with low detectability of modeled SOC changes. 3.3. The Influence of Time Span and Sampling Times Table 5 shows no substantial differences in the mean values of signal, noise, or SNR for δ ^ across different time intervals, regardless of the approach used. Notably, SNR values based on nine‐year SOC observation pairs are not higher than those derived from three‐year intervals. This is also observed in the δ reported by LUCAS ( Supplementary Notebook 04b ). TABLE 5. Summary statistics of signal, noise, and SNR for δ across different time intervals, based on predictions from state‐first and change‐first approach. Category Δ T = 3 Δ T = 6 Δ T = 9 State‐first Change‐first State‐first Change‐first State‐first Change‐first Count 1811 1811 1800 1800 1789 1789 Signal 5.12 2.47 4.97 3.34 5.56 3.23 Noise 63.09 31.74 64.35 32.21 64.30 32.17 SNR 0.09 0.06 0.09 0.08 0.10 0.08 Open in a new tab Compared to δ , which is derived from paired SOC observations between 2009 and 2018, the estimation of β over the same period incorporates an additional intermediate SOC observation. When comparing the SNR of δ ^ and β ^ across this nine‐year span, β ^ (0.11 for the state‐first and 0.16 for the change‐first approach) consistently exhibits slightly higher SNR values than δ ^ (0.10 for the state‐first and 0.08 for the change‐first approach), regardless of the approach used. 3.4. SNR at Larger Spatial Support The correlation function of the standardized errors from c ^ was estimated using variograms to support the calculation of uncertainty in spatial aggregates. Separate variograms were derived for Germany and Spain (Figure 8 ). While both exhibit similar nugget values (variance at zero distance; often attributed to measurement error), the variogram for Germany has a noticeably larger spatial range parameter (about 10.60 km vs. 5.99 km for Spain), indicating that prediction errors in Germany remain spatially correlated over longer distances. FIGURE 8. Open in a new tab Variograms and corresponding correlation functions of standardized prediction errors for c in Germany (left, nugget: 0.72, range: 10.60) and Spain (right, nugget: 0.74, range: 5.99). Figure 9 ‐left panel shows that for δ ^ , both the signal and the noise generally decline as the spatial aggregation support increases, except for a slight signal increase in Spain between 40 km and 200 km. Overall, the reduction in noise is greater than the reduction in signal, resulting in a steady increase in SNR in both countries. The spread of both signal and noise also narrows with increasing aggregation size, whereas the spread in SNR widens. The SNR increase with spatial aggregation is less pronounced in Germany than in Spain. FIGURE 9. Open in a new tab Box plots of signal, noise, and SNR for δ ^ and β ^ in Germany and Spain for spatial aggregation unit sizes of 2.5 km, 10 km, 40 km and 200 km. Boxes represent the interquartile range, and whiskers indicate the remainder of the distribution. Note the difference in the y‐axes. For β ^ , both the signal and noise also decrease as the aggregation size increases. However, the rates of decrease are similar, so the SNR improves only marginally. In Germany, the SNR remains relatively constant across all aggregation levels, while in Spain it increases only slightly. As the aggregation unit size increases, δ ^ values tend to converge toward zero, yet many grid cells in both countries exhibit increasing SNR values (Figure 10 ). The spatial distribution of δ ^ and its SNR from 2009 to 2018 indicates that the SNR increase with larger units is more pronounced in Spain, where several grid cells exceed an SNR of one. In contrast, although SNR also increases in Germany, most grid cells remain below one. FIGURE 10. Open in a new tab Spatial distribution of δ ^ (left two columns) and corresponding SNR (right two columns) from 2009 to 2018 in Germany and Spain, quantified using the state‐first approach. Results are shown across increasing spatial aggregation scales: 1 km pixel‐level predictions, and aggregated units of 2.5 km, 10 km, 40 km and 200 km. Map lines delineate study areas and do not necessarily depict accepted national boundaries. 4. Discussion 4.1. Poor Accuracy in Predicting SOC Change We evaluated two approaches to modeling SOC change— state‐first and change‐first —with δ and β as target variables. Across all four scenarios, prediction accuracy for SOC change remained low (Figure 3 ). Most existing DSM SOC time series adopt the state‐first approach (Gasch et al. 2015 ; Hengl et al. 2017 ; Poggio et al. 2021 ; Heuvelink et al. 2021 ; Li et al. 2022 ; Helfenstein et al. 2024 ; van Wesemael et al. 2024 ; Szatmári, Pásztor, et al. 2024 ; Tian, de Bruin, et al. 2025 ), which are often used for change detection and driver analysis. However, consistent with findings from the Netherlands (Helfenstein et al. 2024 ), the accuracy of SOC change modeling remains poor despite high accuracy in c ^ . A key limitation of the state‐first approach is that it treats predictions for different time instances as independent variables, ignoring temporal autocorrelation in both SOC and environmental covariates (Poggio et al. 2021 ; Szatmári, Pásztor, et al. 2024 ; Tian, de Bruin, et al. 2025 ). As a critical component of SCORPAN (McBratney et al. 2003 ), the time factor is therefore not well integrated into current state‐first space–time DSM frameworks. Efforts to address this gap include incorporating the age of observations as a covariate (Batjes et al. 2020 ), using lagged vegetation indices that weight past conditions to reflect delayed SOC responses (Heuvelink et al. 2021 ), or explicitly incorporating past land use changes in the input covariates (Helfenstein et al. 2024 ). However, constructing such covariates is often impractical at large scales. Moreover, spatial variability in SOC typically outweighs subtle temporal variation, limiting the benefits of such preparation (Poggio et al. 2021 ; Tian, de Bruin, et al. 2025 ). This pattern is also reflected in the exploratory feature contribution analysis (see Supplementary Notebook 14 ), where static features exert a stronger overall influence on model predictions than temporal variables. As a result, models tend to overfit spatial patterns while underrepresenting temporal dynamics, projecting spatial relationships—such as higher SOC under denser vegetation—onto the temporal dimension, even when actual SOC changes occur more slowly than vegetation. The change‐first approach takes a different strategy by modeling SOC change directly. In theory, this offers two advantages: (1) correlated noise in SOC values may cancel out when computing change, and (2) delayed environmental effects may be better captured by directly linking environmental dynamics to SOC change. However, as shown in Figure 3 , SOC change prediction accuracy did not improve as expected. Although the predictive accuracy of β ^ was slightly higher than that of δ ^ , both appeared largely random when validated against repeated measurements, with predictions clustering around zero. This reflects not only the model's limited ability to resolve change, but also the inherent nature of SOC as a slow‐changing property, for which true changes occur over long time scales and are typically small over short intervals (Smith 2004 ; Poeplau and Don 2013 ; Gubler et al. 2019 ). This finding aligns with Broeg et al. ( 2024 ), who attempted to predict β using a much longer (35‐year) cropland SOC time series in Bavaria. While their model struggled to accurately estimate β , it succeeded in classifying SOC change direction (increase, decrease, stable). In contrast, with only the nine‐year time span of LUCAS observations, our study shows that β ^ predictions could not even resolve the direction of change. 4.2. SNR as a Predicted SOC Change Assessment Metric Several studies have applied the concept of SNR to describe challenges in detecting SOC change. For instance, Stevens et al. ( 2008 ), Croft et al. ( 2012 ) reported low SNR in EO reflectance data caused by illumination, terrain, and atmospheric effects, while Paustian et al. ( 2019 ) highlighted that SOC changes are often masked by natural spatial variability. Broeg et al. ( 2024 ) quantified SNR by defining the signal as a temporal variation in measured SOC over a defined time interval and the noise as the corresponding spatial prediction error. Although the specific definitions and calculations of SNR differ across studies, our findings are consistent with those of Broeg et al. ( 2024 ). In both cases, SNR values rarely exceed one over short time spans, surpassing this threshold only after approximately 25 years of observation, and SNR is higher for measurements with higher SOC contents ( > 15 g/kg). This convergence underscores the robustness of SNR as a diagnostic metric and the importance of long‐term monitoring for capturing SOC change. Large‐scale SOC change assessments are often based on prediction time series using the state‐first approach, as repeated measurements are typically unavailable. In such contexts, SNR provides a practical and interpretable means of quantifying the detectability of predicted changes when direct validation is infeasible. It can therefore serve as a readily available indicator of the reliability of SOC change estimates in widely used SOC mapping products. When repeated measurements are available, change‐first approaches become feasible, allowing overall error metrics to be computed by directly comparing predicted and observed changes. However, these metrics remain global summaries and are limited to locations where measurements exist. By contrast, because SNR can be computed for every individual prediction, it can function as a spatially explicit quality indicator across the full prediction domain. In this sense, SNR complements traditional accuracy metrics by providing additional insight into the reliability and interpretability of modeled SOC change, but it should not be viewed as a replacement for predictive accuracy. Accuracy metrics evaluate model performance against observations, whereas SNR serves as an internal diagnostic that indicates whether predicted changes exceed the model's own uncertainty. A key prerequisite for meaningful SNR assessment is valid uncertainty estimation. In this study, all three models— R F c , R F δ , and R F β —produced reasonably calibrated uncertainty intervals, as demonstrated by both coverage (Figure 4 ) and sharpness ( Supplementary Notebook 15 ). Under this condition, SNR provides information on whether the model itself considers a predicted change to be distinguishable from its associated uncertainty, independent of whether the prediction is ultimately accurate. We further compared SNR with predictive accuracy at the point level ( Supplementary Notebook 12 ). Only weak correlations were observed between SNR and relative absolute error (Spearman's ρ = 0.17 ) for both δ ^ and β ^ . At first glance, this weak relationship, together with the generally poor predictive accuracy, may raise questions about the validity of the proposed SNR framework. However, across all modeling strategies, we consistently observe both low predictive accuracy and low SNR. These outcomes are mutually consistent, indicating that SOC change signals are dominated by noise and therefore offer limited interpretive value. As soil monitoring networks such as LUCAS expand, national soil surveys continue, and EO and ancillary datasets improve, longer and denser time series are expected to support more robust models and higher signal‐to‐noise ratios. Under such conditions, SNR may align more closely with external accuracy measures and provide a richer basis for investigating where, why, and under which conditions SOC change becomes more detectable. 4.3. Factors Influencing SNR When comparing SNR between δ ^ and β ^ for the same time span, we observed higher SNR values for β ^ , which incorporates a third SOC observation. In contrast, signal, noise, and SNR of δ remained largely stable despite increasing intervals between paired two observations (Table 5 ). This highlights the importance of not only long‐term monitoring, but also adequate temporal sampling density. Using data from Denmark's long‐term Soil Monitoring Network, Harbo et al. ( 2023 ) showed that short‐term fluctuations can mask long‐term SOC trends, making it risky to estimate δ from only two observations. Similar short‐term fluctuations are also evident in the LUCAS repeated measurements Supplementary Notebook 11b , which can be partially smoothed during trend fitting when multiple time steps are available. Consistently, Broeg et al. ( 2024 ) showed that SNR increases with longer observation periods. The availability of LUCAS 2022 data will allow further assessment of this effect. Compared to the state‐first approach, change‐first models generally produced more variable SNR values across land‐cover series for both δ and β . While both approaches yield variation in predicted changes across land‐cover contexts, the change‐first approach exhibits more stable uncertainty estimates across land‐cover types, suggesting the presence of a baseline noise level in its predictions. As a result, differences in SNR across land‐cover series are more strongly driven by variation in the predicted signal than by changes in uncertainty. For specific land‐cover contexts—such as Woodland for δ (Figure 6 ) and both Woodland and Grassland for β (Figure 7 )—the change‐first approach generally yields higher SNR values than the state‐first approach, whereas lower SNR values are observed for other land‐cover settings. Notably, these patterns appear to be more closely related to the presence of Woodland or Grassland in the land‐cover history than to the occurrence of land‐cover transitions themselves. For example, sites classified as Grassland throughout the entire observation period show higher SNR values than sites experiencing Grassland–Cropland or Cropland–Grassland transitions (Figure 7 ), despite the latter being expected to exhibit larger SOC changes. Although increasing vegetation cover can mask soil signals and reduce detectability from a sensing perspective, dense vegetation is often associated with higher SOC levels, which may counterbalance this effect. This observation is consistent with the findings of Broeg et al. ( 2024 ), who reported that higher SNR is associated with larger background SOC signals, which are commonly found in Woodland and Grassland systems. 4.4. Spatial Aggregation Improves SNR, but Not Always The results confirm that spatial aggregation can improve SNR (Figure 9 ). For δ , noise decreased more rapidly than signal with increasing aggregation size, resulting in a clear increase in SNR. The reduction in noise is primarily a statistical effect: when aggregating over larger spatial units, random and spatially uncorrelated errors tend to cancel out (Smith 2004 ; Wadoux and Heuvelink 2023 , 2025 ). In contrast, the aggregated signal becomes smaller and less variable as it is dominated by the prevalence of small or near‐zero changes within the aggregation unit. For β , both signal and noise decrease with increasing aggregation size at relatively similar rates, leading to largely stable SNR values, with only a slight increase at the largest aggregation scale (200 km 2 ). This difference likely reflects the fact that β is estimated from multiple time steps, which already suppresses part of the short‐term variability and measurement noise during the trend‐fitting stage. As a result, spatial aggregation provides less additional noise reduction for β ^ than for δ ^ . Despite the weaker aggregation effect, β ^ generally exhibits higher SNR values than δ ^ , indicating greater detectability of trend‐based change estimates across spatial scales. Another factor influencing SNR behavior is the spatial domain being mapped—specifically, its environmental conditions, the underlying pattern of pixel‐level SNR, and the spatial structure of model errors. In our study, the improvement was more pronounced in Spain (Figures 9 and 10 ). Although SNR also increased with aggregation in Germany, the improvement was smaller and insufficient to exceed one, even at 200 km 2 support. This contrast primarily reflects environmental and land‐cover differences: Spain's more sparsely vegetated landscapes produced smaller signal magnitudes but also lower uncertainties than Germany (Figures 6 and 7 ). Furthermore, the longer variogram range observed in Germany suggests that spatially correlated errors persist over larger distances, limiting the benefits of aggregation (Wadoux and Heuvelink 2023 ). The choice of aggregation size and unit also influences SNR behavior through the Modifiable Areal Unit Problem (MAUP) (Openshaw and Taylor 1979 ; Dark and Bram 2007 ). The results revealed the scale effect of the MAUP, where analytical outcomes vary with aggregation level. While SNR generally increased with coarser resolution, the trend was not strictly monotonic—for instance, in parts of northeastern Spain, SNR rose at 40 km 2 but declined again at 200 km 2 as low‐SNR areas became dominant within the larger grids. The zonal effect of MAUP—where different zoning schemes (e.g., aggregation by regular grids or administrative units) yield distinct spatial patterns—also needs consideration, although it was not examined in this study. When pixel‐level SNR is homogeneous within a unit, aggregation reduces noise and preserves the signal, thereby improving SNR, as observed in parts of northwestern and central Spain. In contrast, in heterogeneous areas, strong local signals may be masked by surrounding low‐SNR regions, limiting the benefits of aggregation. Aggregation strategies should also consider their practical relevance for MRV, which depends on contextual factors such as land use, management practices, and local regulations that vary spatially (Smith 2008 ; Gocke et al. 2023 ). For instance, management may differ considerably between neighboring farms, yet in some areas, practices are relatively homogeneous due to shared traditions or policy requirements. In this study, we used regular grids for demonstration purposes, but in practice, aggregation units should reflect local land‐use and management contexts. For example, in the Netherlands, farmers are legally required to grow cover crops after maize cultivation (Fan et al. 2020 ), resulting in relatively uniform management across large areas. In such cases, aggregation at administrative or farm levels may be more meaningful than uniform grid cells. However, this improvement comes with a trade‐off: increasing aggregation enhances detectability but reduces spatial detail and local interpretability. Although higher SNR with aggregation is statistically expected, it does not necessarily imply improved model skill. Depending on the aggregation scheme, aggregation may yield little benefit or lead to apparent improvements less meaningful or representative. Ultimately, incorporating domain knowledge, local context, and intended use into aggregation design is necessary to ensure the interpretability and policy relevance of SOC change assessments. 4.5. Strengths and Challenges of Using LUCAS for SOC Change Modeling LUCAS forms the cornerstone of this study. As of this writing, three survey rounds (2009/2012, 2015, and 2018) are publicly available, covering a total of 9 years—still a relatively short period for reliable modeling of SOC change using ML and EO data alone. In Denmark, decadal sampling has been suggested as sufficient to capture long‐term SOC trends, as short‐term fluctuations may not reflect meaningful change (Harbo et al. 2023 ). The upcoming release of LUCAS 2022 will extend the record beyond a decade and provide an additional time step for trend analysis, increasing the potential for detecting long‐term SOC dynamics. While LUCAS is an exceptional resource, several inherent characteristics limit its suitability for SOC change modeling. Designed to represent all of Europe with roughly 20,000 sites (Orgiazzi et al. 2018 ), it effectively captures broad spatial patterns. However, this continental scope can mask subtle local or temporal variations in SOC, particularly in large‐scale mapping frameworks where a single generalized model is applied across all locations and time steps (Lemercier et al. 2022 ). In addition, the relatively low sampling density of such continental surveys restricts their suitability for fine‐scale spatial aggregation, as fitting spatial variograms requires sufficient site density at short distances. For this reason, denser national surveys were used in this study to derive standardized error variograms for the spatial aggregation analysis. Coordinating sampling across many member states and survey years also poses challenges for maintaining consistency. Although all samples are analyzed using standardized methods, variations in both laboratory and field implementations may still introduce non‐negligible differences—particularly given the subtle SOC changes being monitored. Laboratory‐related inconsistencies, for example, have been reported in alkaline soils, where higher analytical uncertainties have been observed in SOC measurements (Schneider et al. 2021 ). Field‐related inconsistencies are also reflected in the extreme apparent SOC changes observed (Figure 3 ). In some cases, δ reaches 400 g/kg and β up to 50 g / kg / year —values that far exceed realistic expectations (Poeplau and Don 2013 ; Gubler et al. 2019 ). Many of these outliers occur in Woodland areas, where inconsistent separation of organic litter from mineral soil layers is suspected (Hiederer 2020 ; Ziche et al. 2022 ). Another likely source of such anomalies is the destructive nature of soil sampling: exact resampling of the same point is impossible. Although LUCAS mitigates this through composite sampling, the official data evaluation report (Hiederer 2020 ) emphasizes that repeated samples should be interpreted as representative of a small area rather than a fixed point. Minor spatial displacement between survey rounds can therefore generate apparent SOC differences unrelated to temporal change, particularly where spatial heterogeneity exceeds temporal variability (Poeplau et al. 2022 ; Broeg et al. 2024 ). Another source of uncertainty lies in the sampling method itself. In LUCAS, soils are sampled with a spade. While this approach facilitates continent‐wide soil surveys at relatively low cost, it raises concerns about whether spade‐based samples accurately represent the upper 0–20 cm layer in certain land uses. Comparative studies from Switzerland indicate that hand‐auger sampling provides greater precision, particularly in Grassland and forest soils with pronounced SOC depth variation (Fernández‐Ugalde et al. 2020 ). Furthermore, repeated LUCAS soil measurements are currently limited to the 0–20 cm topsoil layer. Although this is where most SOC change typically occurs, management effects in croplands have been reported to extend to depths of up to 50 cm (Skadell et al. 2023 ). Several refinements could further enhance LUCAS's capacity for SOC change modeling and monitoring. Permanent site marking or more detailed and standardized location referencing could help reduce relocation uncertainty. Likewise, replacing spade‐based sampling with hand augers would improve depth control and the representativeness of samples within the 0–20 cm layer. LUCAS already provides detailed survey‐level quality reports—an excellent practice that promotes transparency. Extending this to the individual record level, for example by including per‐sample quality flags (analogous to pixel‐level quality indicators in remote‐sensing datasets), would further enhance data consistency and interoperability. National soil surveys—such as those from Spain and Germany used here—complement LUCAS by adding spatial detail and sampling depth, and stronger coordination with national monitoring networks and EO datasets could further improve temporal consistency and model robustness (Froger et al. 2024 ). While some refinements could increase precision, they would also raise costs. Therefore, we recommend prioritizing denser and more frequent resampling of existing LUCAS sites to build consistent long‐term time series, supported by improved and standardized metadata for uncertainty assessment. Taken together, these recommendations are not criticisms of LUCAS but reflect the challenges of applying a continental‐scale survey to high‐resolution SOC change modeling. LUCAS was designed to provide harmonized, long‐term environmental data rather than detailed modeling inputs and already achieves a strong balance between spatial coverage, cost, and feasibility. The program has continuously evolved based on lessons from previous rounds, and the recent inclusion of subsurface layers (e.g., 20–30 cm) marks a valuable step toward capturing vertical SOC trends. Its continuation remains vital for SOC dynamics modeling and will further strengthen the empirical basis for model calibration and improve SNR. 4.6. Implications for SOC Monitoring and DSM Applications Our findings support the conclusions of Broeg et al. ( 2024 ): reliable SOC change detection at fine spatial resolutions (e.g., point, site, or pixel level) remains difficult when relying on EO data and ML methods alone. High prediction uncertainty, limited time‐series length, and the inherently slow dynamics of SOC all contribute to this limitation. Making these constraints explicit and measurable is a necessary step toward more transparent interpretation of DSM‐based SOC change products. In this context, SNR provides a quantitative diagnostic by explicitly relating predicted change magnitude to its associated uncertainty. Low‐SNR values should therefore be interpreted as a warning against over‐interpretation of modeled SOC change, indicating that predictions are likely dominated by noise. SNR can nevertheless be used in spatial aggregation analyses, where aggregation may improve detectability. In addition, spatial patterns of SNR help identify regions with relatively higher or lower detectability, and locations or aggregation units with higher SNR may serve as starting points for investigating environmental or land‐use conditions associated with more detectable SOC change signals. Improving SOC change detectability will depend heavily on the quality of data. To achieve higher SNR, longer, consistent observational datasets—like those from the LUCAS program—are needed to improve the training, validation, and understanding of SOC change models. Beyond extending temporal depth or increasing spatial support, additional strategies can improve SNR, including the use of more relevant, higher‐quality covariates and advances in model architectures (Heuvelink et al. 2021 ; Tian, de Bruin, et al. 2025 ). Hybrid approaches that combine the strengths of ML models and domain‐specific process knowledge have emerged as a promising direction (Reichstein et al. 2019 ). Notable examples include constraining numerical models with domain knowledge (De Rosa et al. 2024 ), coupling process‐based and ML models through loss functions (Zhang, Heuvelink, et al. 2024 ), replacing sub‐modules of complex process‐based models with ML components (Liu et al. 2024 ), and augmenting limited temporal SOC datasets with synthetic data generated by process‐based models (Zhang, Heuvelink, et al. 2024 ). That said, time series of high‐resolution DSM products remain valuable. In fact, c ^ often yield SNR larger than one, demonstrating their value in capturing spatial SOC patterns. Moreover, this SNR evaluation framework can be adapted to work with typical DSM outputs. These predictions can also be aggregated to coarser spatial units, where the resulting SOC changes tend to exhibit higher SNR, indicating more certain and interpretable trends. Given the limited availability of repeated SOC measurements, especially in regions without long‐term monitoring programs like LUCAS, the state‐first approach will likely remain the most viable strategy for detecting and reporting SOC change in the near future. As long as corresponding SNR values are reported alongside SOC change estimates, this approach still support reliable interpretation of model outputs and help identify areas of meaningful change. We advocate for more inclusion of SNR metrics in SOC change estimate efforts to better account for uncertainty and guide data‐informed decisions. 4.7. Limitations The SNR analysis framework developed in this study evaluates the reliability of SOC dynamics modeling by comparing predicted SOC change magnitudes (signal) against model‐derived uncertainty (noise). Since both components are produced by the ML model, the validity of the resulting SNR values depends on the accuracy of both the predictions and the associated uncertainty estimates. In this study, we employed RF and QRF, which are commonly used in ML‐based DSM to generate predictions with uncertainty (Wadoux et al. 2020 ). While we expect our conclusions to hold for other models with robust uncertainty quantification, further work is needed to confirm their generalizability. Additionally, while we included the trend‐based metric β in our analysis, it was estimated from only three SOC measurements over time. Such limited temporal depth introduces uncertainty in trend estimation. Given current data availability, this remains the most feasible approach. However, the forthcoming LUCAS 2022 dataset is expected to provide an additional time step, enabling more reliable trend estimates and reinforcing the findings of this study. Although predicted changes inherently contain both magnitude and direction, the SNR metric, as implemented in this study, quantifies magnitude‐based detectability. With respect to direction, under the assumption of symmetric prediction uncertainty, a sufficient but not necessary condition for directional inference is an SNR value exceeding 0.5, which implies that the prediction distribution lies entirely on one side of zero change. In this study, however, SNR values remain below this level for the vast majority of predictions, substantially limiting the practical relevance of directional inference under current data availability. When the conditional prediction distribution is asymmetric, directionality cannot be inferred from a single summary measure of uncertainty; upper and lower uncertainty bounds must be considered separately. Under such conditions, change direction may still be resolvable even when overall SNR remains low. While the direction of change is an important aspect of SOC dynamics, its reliable assessment requires a dedicated methodological framework, potentially involving classification‐based models and explicit treatment of directional uncertainty, which lies beyond the scope of the present study and could be addressed in future work. While this study provides a quantitative demonstration of how spatial aggregation can improve the SNR, the analysis was intentionally simplified to illustrate the concept rather than to exhaustively evaluate all possible aggregation scenarios. The results should therefore be viewed as a demonstration of the general potential to identify SOC trends more reliably at larger spatial scales, rather than as an attempt to quantify national or regional SOC changes (e.g., in Germany or Spain). Land use and management—key drivers of SOC dynamics—were not explicitly considered in this analysis due to the lack of harmonized, temporally resolved management data at the continental scale. Land‐use data at large scales, although increasingly available through ML‐based EO products, remain subject to uncertainty—just like to that of SOC maps derived from DSM. Data on management practices are even scarcer: existing datasets are typically static in time (Sandström et al. 2023 ) or available only at coarse spatial resolution (Lampach et al. 2025 ). The omission of these factors means that part of the observed noise, particularly at the point scale, may reflect unaccounted management heterogeneity. Moreover, the spatial aggregation analysis was restricted to the state‐first approach using the R F c model. Extending the framework to R F δ or R F β models would in principle be possible, but would require denser temporal datasets to support reliable variogram fitting—data that are not yet available. Future work should aim to integrate contextual information—such as land use and management practices—into both SOC dynamic modeling and spatial aggregation analyses. Incorporating these factors would strengthen the mechanistic understanding of SOC change and enhance the policy relevance of results. Further research should also systematically evaluate how different aggregation schemes—varying in both scale and zoning design—affect SNR behavior across contexts (e.g., farm management, policy reporting, and global monitoring). Such analyses would help clarify the mechanisms by which aggregation enhances or diminishes detectability and provide empirical guidance for designing effective MRV systems. Without detailed land‐use and management information or a systematic assessment of aggregation effects, the interpretability of aggregated SNR values—and thus their usefulness for understanding underlying SOC processes—remains limited. 5. Conclusions With SNR consistently below one and low prediction accuracy, our study demonstrates that, under current data availability, modeling SOC dynamics at the site level using ML and EO data alone remains infeasible. The limited number and frequency of repeated SOC observations constrain the ability to capture meaningful temporal trends. Both the state‐first and change‐first approaches exhibit similarly poor predictive performance, although the change‐first method provides more consistent yet typically narrower uncertainty estimates. While a longer monitoring period does not necessarily increase SNR, incorporating more time steps into SOC change quantification tends to improve it. Higher SNR values are also observed in land covers with greater vegetation content, such as Woodland and Grassland . SNR, as an internal metric for assessing the detectability of modeled SOC changes, cannot replace independent validation. It offers a valuable indicator of change detectability, given the scarcity of repeated SOC observations and can be applied to evaluate the prevailing practice of mapping SOC dynamics through the state‐first approach. When repeated SOC measurements become available, SNR could complement traditional performance metrics by providing spatially explicit insights into the reliability of modeled SOC changes. We argue that SNR should be routinely reported when interpreting or claiming SOC change from model predictions. Without explicitly accounting for uncertainty, such maps risk overstating their reliability for MRV purposes. Finally, our analysis confirms that spatial aggregation improves SNR, with the degree of improvement depending on both the spatial structure of model predictions and the choice of aggregation units. This suggests that, while site‐level assessments remain challenging, reliable SOC change modeling is still achievable at larger scales. This offers a promising pathway for regional‐scale monitoring and reporting. Nonetheless, careful selection of aggregation scales and units is required to ensure that reported SNR values remain meaningful and representative. Author Contributions Xuemeng Tian: conceptualization, methodology, software, data curation, validation, writing – review and editing, visualization, investigation, writing – original draft, formal analysis. Sytze de Bruin: conceptualization, writing – review and editing, supervision, investigation, visualization, formal analysis, methodology. Florian Schneider: conceptualization, writing – review and editing, writing – original draft, supervision, investigation, data curation. Martin Herold: conceptualization, writing – review and editing, supervision. Kirsten de Beurs: conceptualization, writing – review and editing, supervision. Conflicts of Interest The authors declare no conflicts of interest. Acknowledgments This research was funded by the European Union under the Horizon Europe programme through AI4SoilHealth (Grant Agreement No. 101086179) and Intergenerational Open Geospatial Carbon Registry (Grant Agreement No. 101218854). The authors acknowledge the use of ChatGPT (OpenAI) to assist in improving the clarity and readability of the manuscript. The authors are grateful to the Thünen‐Institut für AK Agrarklimaschutz for providing access to the German soil laboratory data Bodenzustandserhebung Landwirtschaft (BZE‐LW). Within the Thünen‐Institut für AK Agrarklimaschutz, the first author extends sincere thanks to all colleagues at the Institute of Climate‐Smart Agriculture, as well as to Tom Broeg from the Institute of Farm Economics, for open, constructive, and critical discussions during a research visit, which greatly helped shape this study. Their academic support and warm hospitality are deeply appreciated. The first author also wishes to thank Professor David Rossiter for his generous support, insightful discussions on SOC change detection using DSM products, and contributions to revising and refining the manuscript. Data Availability Statement The data and code supporting the findings of this study are openly available via Zenodo at https://doi.org/10.5281/zenodo.18706448 . Soil data from Germany (BZE‐LW) are openly available at https://doi.org/10.3220/DATA20200203151139 . Soil data from Spain (Parcelas COS and Parcelas INES) are openly available at https://www.miteco.gob.es/content/dam/miteco/es/biodiversidad/servicios/banco‐datos‐naturaleza/2‐cos/bbdd‐cos.zip . LUCAS Soil Survey data are available from the European Soil Data Centre (ESDAC) for the respective survey years: 2009 and 2012 at https://esdac.jrc.ec.europa.eu/content/lucas‐2009‐topsoil‐data , 2015 at https://esdac.jrc.ec.europa.eu/content/lucas2015‐topsoil‐data , and 2018 at https://esdac.jrc.ec.europa.eu/content/lucas‐2018‐topsoil‐data . References Araza, A. , de Bruin S., Herold M., et al. 2022. “A Comprehensive Framework for Assessing the Accuracy and Uncertainty of Global Above‐Ground Biomass Maps.” Remote Sensing of Environment 272: 112917. [ Google Scholar ] Baghdadi, N. , Cerdan O., Zribi M., et al. 2008. “Operational Performance of Current Synthetic Aperture Radar Sensors in Mapping Soil Surface Characteristics in Agricultural Environments: Application to Hydrological and Erosion Modelling.” Hydrological Processes 22, no. 1: 9–20. [ Google Scholar ] Batjes, N. H. , Ribeiro E., and Van Oostrum A.. 2020. “Standardised Soil Profile Data to Support Global Mapping and Modelling (Wosis Snapshot 2019).” Earth System Science Data 12, no. 1: 299–320. [ Google Scholar ] Bauer‐Marschallinger, B. , Freeman V., Cao S., et al. 2018. “Toward Global Soil Moisture Monitoring With Sentinel‐1: Harnessing Assets and Overcoming Obstacles.” IEEE Transactions on Geoscience and Remote Sensing 57, no. 1: 520–539. [ Google Scholar ] Bellamy, P. H. , Loveland P. J., Bradley R. I., Lark R. M., and Kirk G. J.. 2005. “Carbon Losses From All Soils Across England and Wales 1978–2003.” Nature 437, no. 7056: 245–248. [ DOI ] [ PubMed ] [ Google Scholar ] Broeg, T. , Don A., Wiesmeier M., Scholten T., and Erasmi S.. 2024. “Spatiotemporal Monitoring of Cropland Soil Organic Carbon Changes From Space.” Global Change Biology 30, no. 12: e17608. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Chen, D. , Chang N., Xiao J., Zhou Q., and Wu W.. 2019. “Mapping Dynamics of Soil Organic Matter in Croplands With Modis Data and Machine Learning Algorithms.” Science of the Total Environment 669: 844–855. [ DOI ] [ PubMed ] [ Google Scholar ] Croft, H. , Kuhn N., and Anderson K.. 2012. “On the Use of Remote Sensing Techniques for Monitoring Spatio‐Temporal Soil Organic Carbon Dynamics in Agricultural Systems.” Catena 94: 64–74. [ Google Scholar ] Dark, S. J. , and Bram D.. 2007. “The Modifiable Areal Unit Problem (Maup) in Physical Geography.” Progress in Physical Geography 31, no. 5: 471–479. [ Google Scholar ] De Rosa, D. , Ballabio C., Lugato E., Fasiolo M., Jones A., and Panagos P.. 2024. “Soil Organic Carbon Stocks in European Croplands and Grasslands: How Much Have We Lost in the Past Decade?” Global Change Biology 30, no. 1: e16992. [ DOI ] [ PubMed ] [ Google Scholar ] Don, A. , Seidel F., Leifeld J., et al. 2024. “Carbon Sequestration in Soils and Climate Change Mitigation—Definitions and Pitfalls.” Global Change Biology 30, no. 1: e16983. [ DOI ] [ PubMed ] [ Google Scholar ] Even, R. J. , Machmuller M. B., Lavallee J. M., Zelikova T. J., and Cotrufo M. F.. 2025. “Large Errors in Soil Carbon Measurements Attributed to Inconsistent Sample Processing.” Soil 11, no. 1: 17–34. [ Google Scholar ] Fan, X. , Vrieling A., Muller B., and Nelson A.. 2020. “Winter Cover Crops in Dutch Maize Fields: Variability in Quality and Its Drivers Assessed From Multi‐Temporal Sentinel‐2 Imagery.” International Journal of Applied Earth Observation and Geoinformation 91: 102139. [ Google Scholar ] Fernández‐Ugalde, O. , Jones A., and Meuli R. G.. 2020. “Comparison of Sampling With a Spade and Gouge Auger for Topsoil Monitoring at the Continental Scale.” European Journal of Soil Science 71, no. 2: 137–150. [ Google Scholar ] Fritsch, F. N. , and Butland J.. 1984. “A Method for Constructing Local Monotone Piecewise Cubic Interpolants.” SIAM Journal on Scientific and Statistical Computing 5, no. 2: 300–304. [ Google Scholar ] Froger, C. , Tondini E., Arrouays D., et al. 2024. “Comparing Lucas Soil and National Systems: Towards a Harmonized European Soil Monitoring Network.” Geoderma 449: 117027. [ Google Scholar ] Gasch, C. K. , Hengl T., Gräler B., Meyer H., Magney T. S., and Brown D. J.. 2015. “Spatio‐Temporal Interpolation of Soil Water, Temperature, and Electrical Conductivity in 3d+ t: The Cook Agronomy Farm Data Set.” Spatial Statistics 14: 70–90. [ Google Scholar ] Gocke, M. I. , Guigue J., Bauke S. L., et al. 2023. “Interactive Effects of Agricultural Management on Soil Organic Carbon Accrual: A Synthesis of Long‐Term Field Experiments in Germany.” Geoderma 438: 116616. [ Google Scholar ] Gubler, A. , Wächter D., Schwab P., Müller M., and Keller A.. 2019. “Twenty‐Five Years of Observations of Soil Organic Carbon in Swiss Croplands Showing Stability Overall but With Some Divergent Trends.” Environmental Monitoring and Assessment 191: 1–17. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Harbo, L. S. , Lama R., Lemming C., and Elsgaard L.. 2025. “Importance of Sampling Frequency for the Observed Dynamics of Soc Content in the Danish Long‐Term Monitoring Network.” Geoderma Regional 40: e00931. [ Google Scholar ] Harbo, L. S. , Olesen J. E., Lemming C., Christensen B. T., and Elsgaard L.. 2023. “Limitations of Farm Management Data in Analyses of Decadal Changes in Soc Stocks in the Danish Soil‐Monitoring Network.” European Journal of Soil Science 74, no. 3: e13379. [ Google Scholar ] Harper, K. L. , Lamarche C., Hartley A., et al. 2023. “A 29‐Year Time Series of Annual 300 m Resolution Plant‐Functional‐Type Maps for Climate Models.” Earth System Science Data 15, no. 3: 1465–1499. [ Google Scholar ] Helfenstein, A. , Mulder V. L., Heuvelink G. B., and ten Hack‐Broeke M. J.. 2024. “Three‐Dimensional Space and Time Mapping Reveals Soil Organic Matter Decreases Across Anthropogenic Landscapes in The Netherlands.” Communications Earth & Environment 5, no. 1: 130. [ Google Scholar ] Hengl, T. , de Mens Jesus J., Heuvelink G. B., et al. 2017. “Soilgrids250m: Global Gridded Soil Information Based on Machine Learning.” PLoS One 12, no. 2: e0169748. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Hengl, T. , Nussbaum M., Wright M. N., Heuvelink G. B., and Gräler B.. 2018. “Random Forest as a Generic Framework for Predictive Modeling of Spatial and Spatio‐Temporal Variables.” PeerJ 6: e5518. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Heuvelink, G. B.

2014. “Uncertainty Quantification of Globalsoilmap Products.” GlobalSoilMap: Basis of the Global Spatial Soil Information System: 335–340. 10.1007/978-3-319-63439-5_14. [ DOI ] [ Google Scholar ] Heuvelink, G. B. , Angelini M. E., Poggio L., et al. 2021. “Machine Learning in Space and Time for Modelling Soil Organic Carbon Change.” European Journal of Soil Science 72, no. 4: 1607–1623. [ Google Scholar ] Heuvelink, G. B. M.

2018. Uncertainty and Uncertainty Propagation in Soil Mapping and Modelling, 439–461. Springer International Publishing. [ Google Scholar ] Hiederer, R.

2020. “Data Evaluation of Lucas Soil Survey Laboratory 2009 to 2015 Data.” EUR Report, Publications Office of the European Union, EUR 30092 EN. JRC119881. Ho, Y.‐F. , Grohmann C. H., Lindsay J., et al. 2025. “Global Ensemble Digital Terrain Modeling and Parametrization at 30 m Resolution (gedtm30): A Data Fusion Approach Based on Icesat‐2, Gedi and Multisource Data.” ResearchSquare 13: e19673. 10.7717/peerj.19673. [ DOI ] [ Google Scholar ] Isik, S. , Minarik R., and Hengl T.. 2024. “Continental Europe Surface Lithology Based on Egdi/Onegeology Map at 1:1 m Scale.” Karger, D. N. , Conrad O., Böhner J., et al. 2017. “Climatologies at High Resolution for the Earth's Land Surface Areas.” Scientific Data 4, no. 1: 1–20. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Karger, D. N. , Wilson A. M., Mahony C., Zimmermann N. E., and Jetz W.. 2021. “Global Daily 1 Km Land Surface Precipitation Based on Cloud Cover‐Informed Downscaling.” Scientific Data 8, no. 1: 307. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Kasraei, B. , Heung B., Saurette D. D., Schmidt M. G., Bulmer C. E., and Bethel W.. 2021. “Quantile Regression as a Generic Approach for Estimating Uncertainty of Digital Soil Maps Produced From Machine‐Learning.” Environmental Modelling & Software 144: 105139. [ Google Scholar ] Lal, R.

2004. “Soil Carbon Sequestration to Mitigate Climate Change.” Geoderma 123, no. 1–2: 1–22. [ Google Scholar ] Lampach, N. , Skoien J. O., Ramos H., et al. 2025. “Statistical Atlas of European Agriculture: Gridded Data From the Agricultural Census 2020 and the Spatial Distribution of Cap Contextual Indicators.” Earth System Science Data Discussions 2025: 1–35. [ Google Scholar ] Lehmann, J. , Bossio D. A., Kögel‐Knabner I., and Rillig M. C.. 2020. “The Concept and Future Prospects of Soil Health.” Nature Reviews Earth & Environment 1, no. 10: 544–553. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Lemercier, B. , Lagacherie P., Amelin J., et al. 2022. “Multiscale Evaluations of Global, National and Regional Digital Soil Mapping Products in France.” Geoderma 425: 116052. [ Google Scholar ] Lettens, S. , Van Orshoven J., Van Wesemael B., Muys B., and Perrin D.. 2005. “Soil Organic Carbon Changes in Landscape Units of Belgium Between 1960 and 2000 With Reference to 1990.” Global Change Biology 11, no. 12: 2128–2140. [ DOI ] [ PubMed ] [ Google Scholar ] Li, H. , Wu Y., Liu S., et al. 2022. “Decipher Soil Organic Carbon Dynamics and Driving Forces Across China Using Machine Learning.” Global Change Biology 28, no. 10: 3394–3410. [ DOI ] [ PubMed ] [ Google Scholar ] Liu, L. , Zhou W., Guan K., et al. 2024. “Knowledge‐Guided Machine Learning Can Improve Carbon Cycle Quantification in Agroecosystems.” Nature Communications 15, no. 1: 357. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Lyapustin, A. , and Wang Y.. 2018. “MCD19A2 MODIS/Terra+Aqua Land Aerosol Optical Depth Daily L2G Global 1 km SIN Grid V006.” NASA EOSDIS Land Processes Distributed Active Archive Center. McBratney, A. B. , Santos M. M., and Minasny B.. 2003. “On Digital Soil Mapping.” Geoderma 117, no. 1–2: 3–52. [ Google Scholar ] Meinshausen, N.

2006. “Quantile Regression Forests.” Journal of Machine Learning Research 7, no. 6: 983–999. [ Google Scholar ] Meng, X. , Bao Y., Luo C., Zhang X., and Liu H.. 2024. “Soc Content of Global Mollisols at a 30 m Spatial Resolution From 1984 to 2021 Generated by the Novel ml‐Cnn Prediction Model.” Remote Sensing of Environment 300: 113911. [ Google Scholar ] Openshaw, S. , and Taylor P.. 1979. “A Million or so Correlated Coefficients: Three Experiment on the Modifiable Areal Unit Problem.” Statistical applications in the spatial sciences. Orgiazzi, A. , Ballabio C., Panagos P., Jones A., and Fernández‐Ugalde O.. 2018. “Lucas Soil, the Largest Expandable Soil Dataset for Europe: A Review.” European Journal of Soil Science 69, no. 1: 140–153. [ Google Scholar ] Paustian, K. , Collier S., Baldock J., et al. 2019. “Quantifying Carbon for Agricultural Soil Management: From the Current Status Toward a Global Soil Information System.” Carbon Management 10, no. 6: 567–587. [ Google Scholar ] Poeplau, C. , and Don A.. 2013. “Sensitivity of Soil Organic Carbon Stocks and Fractions to Different Land‐Use Changes Across Europe.” Geoderma 192: 189–201. [ Google Scholar ] Poeplau, C. , Don A., Flessa H., Heidkamp A., Jacobs A., and Prietz R.. 2020. First German Agricultural Soil Inventory – Core Dataset. OpenAgrar. [ Google Scholar ] Poeplau, C. , Prietz R., and Don A.. 2022. “Plot‐Scale Variability of Organic Carbon in Temperate Agricultural Soils—Implications for Soil Monitoring#.” Journal of Plant Nutrition and Soil Science 185, no. 3: 403–416. [ Google Scholar ] Poggio, L. , De Sousa L. M., Batjes N. H., et al. 2021. “Soilgrids 2.0: Producing Soil Information for the Globe With Quantified Spatial Uncertainty.” Soil 7, no. 1: 217–240. [ Google Scholar ] Potapov, P. , Turubanova S., Hansen M. C., et al. 2022. “Global Maps of Cropland Extent and Change Show Accelerated Cropland Expansion in the Twenty‐First Century.” Nature Food 3, no. 1: 19–28. [ DOI ] [ PubMed ] [ Google Scholar ] Reichstein, M. , Camps‐Valls G., Stevens B., et al. 2019. “Deep Learning and Process Understanding for Data‐Driven Earth System Science.” Nature 566, no. 7743: 195–204. [ DOI ] [ PubMed ] [ Google Scholar ] Rogge, D. , Bauer A., Zeidler J., Mueller A., Esch T., and Heiden U.. 2018. “Building an Exposed Soil Composite Processor (Scmap) for Mapping Spatial and Temporal Characteristics of Soils With Landsat Imagery (1984–2014).” Remote Sensing of Environment 205: 1–17. [ Google Scholar ] Sandström, E. , Namasivayam A., Oostdijk S., Scherpenhuijzen N., Debonne N., and Verburg P.. 2023. Land System Map for Europe. Schmidinger, J. , and Heuvelink G. B.. 2023. “Validation of Uncertainty Predictions in Digital Soil Mapping.” Geoderma 437: 116585. [ Google Scholar ] Schneider, F. , Poeplau C., and Don A.. 2021. “Predicting Ecosystem Responses by Data‐Driven Reciprocal Modelling.” Global Change Biology 27, no. 21: 5670–5679. [ DOI ] [ PubMed ] [ Google Scholar ] Sen, P. K.

1968. “Estimates of the Regression Coefficient Based on Kendall's Tau.” Journal of the American Statistical Association 63, no. 324: 1379–1389. [ Google Scholar ] Serrano, L. R. , Ruiz A. M., Llorente R. B., Donoso J. J. S., and Pérez S. M.. 2022. “Mapa Del Carbono Orgánico Del Suelo en España: Estimación a Partir De Los Datos Del Inventario Nacional de Erosión de Suelos.” Ministerio Para la Transición Ecológica y el Reto Demográfico (MITECO), Madrid, España. Los Contenidos de Esta Publicación Podrán Ser Reutilizados Citando la Fuente y la Fecha, en su Caso, de la Última Actualización. Shimada, M. , and Ohtaki T.. 2010. “Generating Large‐Scale High‐Quality Sar Mosaic Datasets: Application to Palsar Data for Global Monitoring.” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 3, no. 4: 637–656. [ Google Scholar ] Skadell, L. E. , Schneider F., Gocke M. I., et al. 2023. “Twenty Percent of Agricultural Management Effects on Organic Carbon Stocks Occur in Subsoils–Results of Ten Long‐Term Experiments.” Agriculture, Ecosystems & Environment 356: 108619. [ Google Scholar ] Smith, P.

2004. “How Long Before a Change in Soil Organic Carbon Can Be Detected?” Global Change Biology 10, no. 11: 1878–1883. [ Google Scholar ] Smith, P.

2008. “Land Use Change and Soil Organic Carbon Dynamics.” Nutrient Cycling in Agroecosystems 81, no. 2: 169–178. [ Google Scholar ] Smith, P. , Soussana J.‐F., Angers D., et al. 2020. “How to Measure, Report and Verify Soil Carbon Change to Realize the Potential of Soil Carbon Sequestration for Atmospheric Greenhouse Gas Removal.” Global Change Biology 26, no. 1: 219–241. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Soils Revealed . 2025. “Soils Revealed: Global Soil Organic Carbon Mapping Platform.” https://soilsrevealed.org/ . Stevens, A. , van Wesemael B., Bartholomeus H., Rosillon D., Tychon B., and Ben‐Dor E.. 2008. “Laboratory, Field and Airborne Spectroscopy for Monitoring Organic Carbon Content in Agricultural Soils.” Geoderma 144, no. 1–2: 395–404. [ Google Scholar ] Stockmann, U. , Adams M. A., Crawford J. W., et al. 2013. “The Knowns, Known Unknowns and Unknowns of Sequestration of Soil Organic Carbon.” Agriculture, Ecosystems & Environment 164: 80–99. [ Google Scholar ] Sun, Q. , Zhang P., Jiao X., et al. 2023. “A Global Estimate of Monthly Vegetation and Soil Fractions From Spatio‐Temporally Adaptive Spectral Mixture Analysis During 2001–2022.” Earth System Science Data Discussions 2023: 1–30. [ Google Scholar ] Szatmári, G. , Laborczi A., Mészáros J., et al. 2024. “Gridded, Temporally Referenced Spatial Information on Soil Organic Carbon for Hungary.” Scientific Data 11, no. 1: 1–13. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Szatmári, G. , Pásztor L., Takács K., Mészáros J., Benő A., and Laborczi A.. 2024. “Space‐Time Modelling of Soil Organic Carbon Stock Change at Multiple Scales: Case Study From Hungary.” Geoderma 451: 117067. [ Google Scholar ] Theil, H.

1950. “A Rank‐Invariant Method of Linear and Polynomial Regression Analysis.” Indagationes Mathematicae 12, no. 85: 173. [ Google Scholar ] Tian, X. , Consoli D., Witjes M., et al. 2025. “Time Series of Landsat‐Based Bimonthly and Annual Spectral Indices for Continental Europe for 2000–2022.” Earth System Science Data 17, no. 2: 741–772. [ Google Scholar ] Tian, X. , de Bruin S., Simoes R., et al. 2025. “Spatiotemporal Prediction of Soil Organic Carbon Density in Europe (2000–2022) Using Earth Observation and Machine Learning.” PeerJ 13: e19605. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] van Wesemael, B. , Abdelbaki A., Ben‐Dor E., et al. 2024. “A European Soil Organic Carbon Monitoring System Leveraging Sentinel 2 Imagery and the Lucas Soil Data Base.” Geoderma 452: 117113. [ Google Scholar ] Vaysse, K. , and Lagacherie P.. 2017. “Using Quantile Regression Forest to Estimate Uncertainty of Digital Soil Mapping Products.” Geoderma 291: 55–64. [ Google Scholar ] Venter, Z. S. , Hawkins H.‐J., Cramer M. D., and Mills A. J.. 2021. “Mapping Soil Organic Carbon Stocks and Trends With Satellite‐Driven High Resolution Maps Over South Africa.” Science of the Total Environment 771: 145384. [ DOI ] [ PubMed ] [ Google Scholar ] Wadoux, A. M.‐C. , and Heuvelink G. B.. 2023. “Uncertainty of Spatial Averages and Totals of Natural Resource Maps.” Methods in Ecology and Evolution 14, no. 5: 1320–1332. [ Google Scholar ] Wadoux, A. M.‐C. , and Heuvelink G. B.. 2025. “Scientists Yet to Consider Spatial Correlation in Assessing Uncertainty of Spatial Averages and Totals.” International Journal of Applied Earth Observation and Geoinformation 139: 104472. [ Google Scholar ] Wadoux, A. M.‐C. , Minasny B., and McBratney A. B.. 2020. “Machine Learning for Digital Soil Mapping: Applications, Challenges and Suggested Solutions.” Earth‐Science Reviews 210: 103359. [ Google Scholar ] Wagner, W. , Bauer‐Marschallinger B., Navacchi C., et al. 2021. “A Sentinel‐1 Backscatter Datacube for Global Land Monitoring Applications.” Remote Sensing 13, no. 22: 4622. [ Google Scholar ] Wan, Z.

2006. “Modis Land Surface Temperature Products Users' Guide.” Institute for Computational Earth System Science, University of California: Santa Barbara, CA, USA, 805:26. Wang, B. , Gray J. M., Waters C. M., et al. 2022. “Modelling and Mapping Soil Organic Carbon Stocks Under Future Climate Change in South‐Eastern Australia.” Geoderma 405: 115442. [ Google Scholar ] Widyastuti, M. T. , Minasny B., Padarian J., and Maggi F.. 2024. “Peatgrids: Mapping Global Peat Thickness and Carbon Stock via Digital Soil Mapping Approach, Dataset.” Yang, R.‐M. , Liu L.‐A., Zhang X., et al. 2022. “The Effectiveness of Digital Soil Mapping With Temporal Variables in Modeling Soil Organic Carbon Changes.” Geoderma 405: 115407. [ Google Scholar ] Yigini, Y. , and Panagos P.. 2016. “Assessment of Soil Organic Carbon Stocks Under Future Climate and Land Cover Changes in Europe.” Science of the Total Environment 557: 838–850. [ DOI ] [ PubMed ] [ Google Scholar ] Zhang, L. , Heuvelink G. B., Mulder V. L., Chen S., Deng X., and Yang L.. 2024. “Using Process‐Oriented Model Output to Enhance Machine Learning‐Based Soil Organic Carbon Prediction in Space and Time.” Science of the Total Environment 922: 170778. [ DOI ] [ PubMed ] [ Google Scholar ] Zhang, T. , Huang L.‐M., and Yang R.‐M.. 2024. “Evaluation of Digital Soil Mapping Projection in Soil Organic Carbon Change Modeling.” Ecological Informatics 79: 102394. [ Google Scholar ] Zhang, Z. , Ding J., Zhu C., et al. 2023. “Historical and Future Variation of Soil Organic Carbon in China.” Geoderma 436: 116557. [ Google Scholar ] Ziche, D. , Grüneberg E., Riek W., and Wellbrock N.. 2022. “Comparison of the Lucas 2015 Inventory With the Second National Forest Soil Inventory.” Thünen Report. Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Data Availability Statement The data and code supporting the findings of this study are openly available via Zenodo at https://doi.org/10.5281/zenodo.18706448 . Soil data from Germany (BZE‐LW) are openly available at https://doi.org/10.3220/DATA20200203151139 . Soil data from Spain (Parcelas COS and Parcelas INES) are openly available at https://www.miteco.gob.es/content/dam/miteco/es/biodiversidad/servicios/banco‐datos‐naturaleza/2‐cos/bbdd‐cos.zip . LUCAS Soil Survey data are available from the European Soil Data Centre (ESDAC) for the respective survey years: 2009 and 2012 at https://esdac.jrc.ec.europa.eu/content/lucas‐2009‐topsoil‐data , 2015 at https://esdac.jrc.ec.europa.eu/content/lucas2015‐topsoil‐data , and 2018 at https://esdac.jrc.ec.europa.eu/content/lucas‐2018‐topsoil‐data . Articles from Global Change Biology are provided here courtesy of Wiley ACTIONS View on publisher site PDF (5.0 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top