More than 85% of the regional and national SCORPAN-type studies identified in dataset 1 were carried out during the last ten years, with the growth rate in publications increasing rapidly during recent years ( Figure 1 ). This can largely attribute to the rapid development of DSM technologies enabling easy access to uniform or standardized environmental covariates needed for SOC estimation at large spatial scales. The goodness of model fit (R²) for estimations of agricultural SOC varied from 0.02 to 0.86, with an averaged value of 0.45 and a median value of 0.47 for all records from the selected regional and national studies ( Table S1 ). Spearman’s rank order correlation between R² and modeling factors indicates that the goodness of model fit decreased with deeper soil sampling depth (Spearman R = -0.327, p < 0.001), larger spatial extent (Spearman R = -0.472, p < 0.001), and smaller grid-distance (Spearman R = 0.240, p < 0.01). The correlation between R² and the number of SCORPAN predictors was positive but not significant (R = -0.063, p > 0.05) ( Table S3 ).

Soil depth has been considered crucial for SOC quantification, and it is well recognized that both surface and subsurface soil layers must be considered for effective assessment (118–120). Numerous works have called for the use of environmental covariates that can predict the vertical distribution of SOC (121–123). Despite this, only 48.1% of the studies (N = 38) quantified SOC beyond 30 cm. Negative correlation between model fit and soil sampling depth suggests that SCORPAN-type models better quantified SOC in surface than subsurface layers ( Table S2 ). While models that included SOC beyond 60 cm depths typically had poorer fits than those considering surface depths ( Figure 2A ) and half of the studies (57.9%) reporting results for separate soil sampling depths had clearly larger RMSE values associated with SOC estimates in deeper soils, some studies still obtained moderate estimates (R² > 0.5) for deeper soil layers ( Table S1 ). Differences among studies from a larger dataset might make it impossible to identify covariates contributing to model success. Use of 3-D modeling approaches to parameterize continuous depth functions has been proposed as one way to efficiently predict subsoil SOC stocks (43, 53, 66, 94, 124–126), but this method does not always improve the accuracy over that obtained with traditional 2.5D modeling frameworks that rely on predetermined depth increments. Ma et al. (127) found “stepped” depth function artifacts can occur for SOC prediction when depth is used as a covariate in tree-based algorithms. Obtaining finer resolution depth measurements to overcome data paucity issues and incorporating key pedology drivers separated for different depth layers in SCORPAN models would be essential for improving SOC estimates from different depths.

Better model fits were typically observed in SCORPAN-type studies that covered smaller spatial extents, particularly when models were based on surface soil layers (< 30 cm) ( Figure 2B ; Table S4 ). This is consistent with the idea that model fits are improved by reducing both vertical and horizontal variability. Our observation that models built with smaller grid-distance (< 500 m), or greater sampling intensity, did not have greater performance ( Figures 2C, D ; Table S4 ) regardless of sampling depth and is inconsistent with the basic rules for statistical sampling. Studies like the work of Minasny et al. (39) find that model’s statistical power increases with better representation of the sampled population. These kinds of geostatistically based conclusions, which are drawn using data mainly from field or local-scale studies with grid spacing up to 1 km, may not apply to regional or national studies summarized by us that consider grid distances ranging from 15 m to more than 100 km. Even though finer spatial resolution should provide more detailed model prediction (128), the lumped effects of variability associated with larger grid distances may have a smoothing effect that reduces complexity and explain why model fits did not decline with grid size expansion (129). Similarly, increasing the total number of predictors involved in SCORPAN-type model building did not have significant influence on model performance based on the studies considered ( Table S3 ). This explains why the assumption that increased availability of data representing SCORPAN covariates would improve model predictions might fail in the situation where added covariates cause artifacts or model-overfitting that do not mechanistically explain the distribution of SOC (130). To further validate the assumption, future work should test the use of different covariate sets to identify the optimized number of predictors needed to capture SOC variability.

Although simple linear regressions were predominantly used, tree-based regression techniques such as Random Forest, Decision Tree, and Cubist Regression Tree were examined in more than half of studies reviewed (53%, N = 42), and these techniques were more often applied in recent publications. Only 22 studies (28%) interpolated regression residuals (unexplained variation) using kriging to account for spatially auto-correlated errors ( Table S1 ), even though spatial information provided by kriging is known to be of vital importance (131, 132). Even fewer studies used data mining techniques such as Support Vector Regression and Artificial Neural Network (10%, N = 8). These machine learning techniques can potentially characterize complicated and indirect relationships between soil properties and environmental covariates (133, 134). Although the reported R² values had mixed results for different modeling approaches, better goodness of fit was generally obtained from tree-based and machine-learning methods than from other techniques ( Table S1 ). To move beyond using a single modeling technique, the study of Tajik et al. (135) reported that the ensembled learning technique, which utilizes multiple learners based on a number of hypotheses, might be more robust at quantifying SOC.

Even though the need for model validation and uncertainty analysis for SOC quantification is well-recognized (136–140), only 49% (N = 39) of the studies we reviewed validated their results independently. This is a slight improvement over earlier reviews (38, 39). Other studies either used cross-validation (37%; N = 29) or did not provide any type of validation (14%; N = 11). We also found that less than half (N = 36, 46%) of the regional and national scale studies included model uncertainty analysis but is more frequently included in more recent studies ( Figure 1 ). Most often, uncertainty is represented by confidence intervals for model predictions that reflect uncertainty associated with model inputs. While data quality is a well-known source of SCORPAN model uncertainty (38, 88), other sources including data interpolation or rescaling (141, 142) and spatial and temporal mismatches (38, 143, 144) were largely unaddressed.

Of the 79 regional and national-scale studies, the most widely used SCORPAN predictors were R (93.7% of the cases; N = 74), C (87.3%; N = 69) and O (86.1%; N = 68), followed by S (57.0%; N = 45) and P (46.8%; N = 37). Less commonly used predictors included N (10.1%; N = 8) and A (5.1%; N = 4) ( Figure 3A ). Over 80% of the regional and national scale studies used covariates representing 3 (25%, N = 20), 4 (34%, N = 27), or 5 (23%, N = 18) SCORPAN predictors, with SCORP, SCOR, COR, and CORP being the most commonly adopted combinations ( Figure 3B ). A smaller percent of studies adopted covariates from two (13%; N = 10) or more than 5 predictors (5%, N = 4). The fact that studies rarely utilize all seven SCORPAN predictors is probably because too many covariates can cause multicollinearity issues (145) or because of parsimony considerations. This trend may alter in the future as more recent studies tend to use a larger number of predictors due to increased data availability. Priorities should be given to covariates that can best represent processes known to influence SOC variability. Because of this, our following texts tried to explain the physical meanings of the SCORPAN covariates based on their frequency of use.

Our finding that relief (R) was the most used predictor ( Figure 3A ) agrees with a review of soil mapping and scale finding that soil maps at large cartographic scales have relied more on topography than other variables (146). The R predictor is critical for water movement and material accumulation (147) and thus can indirectly affect predictors S, C, and O. For example, R covariates were found to influence or correlate with soil moisture conditions, temperature, precipitation, and vegetation patterns (148–151). Many of the studies we reviewed here used multiple covariates from both the ‘primary’ and ‘secondary’ attributes of the R predictor ( Table 2 ). Primary R attributes usually include elevation, slope, aspect, catchment area, and landform curvatures while secondary attributes (e.g., soil water and erosion indices) are derived from them (152–155). Although some R covariates are relatively independent (e.g., Soliveres et al. (156): elevation and topographic index), secondary R covariates may be strongly correlated (e.g., Wilson et al.’s (2000) stream power index and slope length factors). Therefore, future studies may use techniques to identify and exclude R covariates that are repeated and non-important for SOC quantification (157–159).

Climate (C) and organisms (O) covariates are widely used in large spatial extent studies ( Figure 3A ) because both factors are likely to vary to the degree that can significantly influence agricultural SOC (146). Climate is the main driving factor for SOC at the large scale because of its direct influence on SOC decomposition and indirect impacts on other S (i.e. through soil pH and texture), O (i.e. through C inputs), and P (i.e. through mineralogy) covariates (160, 161). Over 60% (N = 42) of the studies that included C covariates only used precipitation and temperature for SOC quantification. The rest of the studies mostly used covariates including potential evapotranspiration, solar radiation, humidity, and vapor pressure (VPD) ( Table 2 ). Predictor O included covariates that relate to both human activities such as agricultural management and vegetation cover and production. Most utilized O covariates are land cover (68.4%, N = 54) and vegetation indices (53.2%, N = 42). Despite the fact that tillage and fertilizer management significantly impact agricultural SOC (162–164) and that the influence of anthropogenic factors on the vertical and horizontal variability of soil properties should be better accounted for within the SCORPAN framework (165, 166), only 6 studies extracted information about manure or fertilizer application and 2 studies extracted crop or rotation type. None of our reviewed studies reported the use of tillage or residue management practices as model covariates.

Soil (S) and parent material (P) covariates were less commonly used than R, O, and C covariates for the large spatial extent studies ( Figure 3A ). While Grunwald (38) found S covariates were used 84% of the cases to estimate various soil attributes and classes, with SOC being used as an S covariate, we found quantification efforts for SOC valued the utility of S factors less (< 60%). Soil classification (25.3%, N = 20) or taxonomic data (24.1%, N = 19), followed by bare soil reflectance (6.3%, N = 5) and soil erosion rate (5.1%, N = 4) were among the more commonly used S covariates ( Table 2 ). Rasmussen et al. (167) proposed that properties such as exchangeable Ca and Fe-oxyhydroxides are more explanatory of SOC due to their associations with SOM stabilization mechanisms. However, those S factors are less used in large-scale models due to a lack of data products. The most widely used P covariates (gamma-ray spectroscopic measures, rock fragments, and geological indices) were employed in less than a quarter of the reviewed studies and so their utilities are difficult to evaluate and compare with other more commonly adopted covariates. This calls for selection and potential incorporation of pedolologically-relevant SCORPAN covariates that can authentically capture distinctive information about SOC variability. Covariates like the subsolum reference groups associated with different weathering stages proposed by Juilleret et al. (168) are examples of P covariates that could improve modeling of subsurface soils. In addition, established soil-geomorphic associations such as simple and complex catenas (169), soil associations, and other aggregations of landscape scale soil patterns, can be configured into elements with embedded soil covariate properties. An example of this approach is illustrated by Atkinson et al. (170) where geomorphon (a geomorphological phonotype) is used for digital geomorphological mapping. They point out that geomorphon feature relevance for defining landscape structure and terrain spatial heterogeneity must be framed in the context of landscape or terrain detail, soil covariate membership, DEM pixel resolution, and user preference. Similarly, Jafari et al. (171) identified geomorphic surfaces and terrain attributes to be effective at capturing spatial patterns in soils. The prospect of using scale-appropriate soil geomorphic units as components of regional-scale investigations, which incorporate continuous variation of soil properties is a potentially lucrative approach for improving modeling of SOC and other soil attributes at the large spatial scale.

Even though SOC is known to vary spatially, predictor N that explicitly describes location or space received relatively little attention ( Figure 2A ). This may be because spatial information can be easily reflected by other SCORPAN predictors such as C or R covariates (145, 172). Studies that considered N covariates generally included spatial coordinates or proximity to objects in the model ( Table 2 ). Arguments in favor of using independent N covariates are that variables are relatively easy to obtain and can be used as proxies for other more complicated variables, and can account for within-grid heterogeneity, correct spatial autocorrelation of model residuals, and explain model uncertainty and spatial patterns not captured by other environmental covariates (173–176). Bjørn Møller et al. (177) proposed the use of coordinates adjusted to oblique angles to alleviate orthogonal artifacts, but this idea has yet to be tested in regional or national scale SCORPAN studies.

Predictor A was even less commonly used than predictor N in studies contained in our dataset ( Figure 2A ). This is consistent with previous reviews finding limited use of SCORPAN covariates tracking temporal trends (20, 38, 39). We found the majority of the studies (89.9%, N = 71) investigated SOC from just one sampling time or assumed that SOC collected within a period of time from years to decades from the national library were constant. Even when studies sampled SOC multiple times to investigate changes in SOC (46, 73, 79, 96, 97, 102), models were built separately for each year rather than based on use of A covariates likely due to the difficulty to resample the same soil profile. The review of Croft et al. (178) showed that RS data may be promising for modeling temporal SOC changes through the monitoring of soil structural changes, soil erosion, agricultural practices in time, but the accuracy of these covariate data obtained at the large spatial extent need to be tested further.

The majority of the regional and national scale studies (76%) directly utilized existing SOC databases for SCORPAN model calibration ( Table 1 ) while smaller percentages of the studies combined survey data with literature values or additional field sampling campaigns (5%). One key issue is that most databases ( Table S5 ) only reported SOC from a single sampling time. Since spatial differences among soil properties can be confounded with temporal changes when soils are collected over relatively long time periods (e.g., decades), it would be challenging to use a single, one-time sampling SOC database to build temporal SCORPAN models. By harmonizing multiple soil databases or using additional measures designed to address spatial variability (179), enhanced SOC datasets may be incorporated into time-space modeling frameworks for calibration to reflect SOC dynamics (23, 38, 180, 181). However, differences among sampling and testing methods used by databases can confound comparisons so soil inventories must use standardized procedures to accurately quantity change in SOC or other dynamic soil properties. Additional opportunities may arise for future work as networks (e.g., The Soils 2026 and Digital Soil Mapping’ initiatives) are established to provide continuous predictions of soil properties and the associated estimates of uncertainty for the U.S. (30).

The S covariates can be accessed from a large number of gridded or point-based databases ( Table S5 ) which most commonly included soil order and series, and measures or pedo-transfer function-based estimates of soil texture, bulk density, pH, available water capacity, and cation exchange capacity, followed by horizon depth, total N, drainage class, and moisture. The gridded soil datasets have spatial resolutions ranging from approximately 30 m to more than 10 km, with the finer resolution databases being mostly interpolated (e.g., gSSURGO (182) from coarser resolution products or estimated with the remote sensing (RS) technique (e.g., POLARIS (183)). The majority of the soil databases report site-based point observations rather than interpolated results and so are often rescaled by the SCORPAN studies to match the mapping unit and resolution of other covariates.

The C covariates including precipitation, temperature, potential evapotranspiration, solar radiation, and VPD have been covered by a large number of databases ( Table S6 ), with the spatial and temporal resolution varying from 1 to 130 km, and from sub-hourly to monthly, respectively. There are generally two different types of databases for C covariates, one type uses directly measured data from stations such as the APHRODITE (184) and the USHCN (185) databases; the other type of datasets for C covariates are interpolated based on measured and RS-derived data and are more widely used by SCORPAN-type studies. For instance, the PRISM dataset that is commonly used by U.S. studies employed elevation for data interpolation (148) and the NASA AIRS dataset relies on the detection of infrared energy emitted from Earth’s surface to derive temperature and water vapor measures. The C covariates are commonly reported at coarser spatial but finer temporal resolution than O or R covariates ( Table S6 ). Unfortunately, C covariates were generally averaged in time as inputs for the SCORPAN models especially when temporal SOC changes were not considered. Ideally, studies could also use variance associated with C covariates as model inputs to couple spatial SOC estimates with temporal drivers.

The O covariates including Normalized Difference Vegetation Index (NDVI), land use or land/vegetation cover (LULC), and net primary production (NPP) can be easily extracted from RS-derived datasets ( Table S6 ). While regional and national scale efforts typically process RS data and derive LULC products along with more detailed crop types (e.g., CDL dataset (186)) annually or every few years, NDVI and other vegetation indices can be estimated at much finer temporal resolutions. Even though only 10% (N = 8) of the reviewed studies utilized NPP as a SCORPAN covariate (Table 2 ), Running et al. (187) predicted that NPP estimates can be reported as frequently as weather products in the future. In contrast, current large-scale datasets addressing O covariates cannot precisely identify where component practices (e.g., fertilizer and manure application (188); tillage and residue management (189)) are applied within agricultural landscapes. Future SCOPRAN-type efforts on SOC modeling should make use of more spatially-resolved management datasets provided that there is a close association between changes in SOC stocks and management practices (190, 191).

The R covariates can be extracted from Digital Elevation Models (DEMs) (192) that are built at the large spatial scales ( Table S6 ). The SCORPAN studies we reviewed most commonly extracted DEMs from the Shuttle Radar Topographic Mission (SRTM) (193), which produced a near-global database with relatively fine spatial resolution (30 m). More recent efforts have applied the light Detection and Ranging (LiDAR) technology to further enhance the data quality and resolution of R covariates by capturing fine terrain features from complicated ecosystems (194, 195). For example, the LiDAR-derived DEMs available for the U.S., U.K., Australia, and Denmark have the finest resolutions among all the R data sources ( Table S6 ).

Although P covariates were not frequently used by studies reviewed, some P covariates such as soil parent material, lithology, bedrock depth, and gamma-ray spectroscopic measures have been included in global or national datasets (Table S6). Likewise, N factors either in the form of coordinates or distance to specified objects is not widely adopted but should be relatively easy to extract from GPS or legacy maps. In contrast, information about A covariates is generally lacking, particularly for legacy data that underpins most soil surveys. Obtaining N factors may need to rely on genesis type (loess, fluvial, aeolian, arid, glacial, and periglacial landforms) or topographic proxies (196).

The advancements of RS have addressed data paucity issues by providing a large number of SOC and SCORPAN covariate datasets ( Tables S5, S6 ). However, not only the spatial and temporal coverages and resolution of the datasets but also their data quality are major concerns for building empirical soil models (61, 197, 198). This calls for future SCORPAN-type efforts to thoroughly compare and check the consistency of the RS datasets before choosing proper sources for SCORPAN covariates. Incorporating temporal RS data as model covariates (e.g., vegetation cover and type, climate factors) can potentially reduce the uncertainty of empirically estimated large spatial scale soil properties (199, 200). Future studies may also consider the use of temporal SOC measures enhanced with RS images (e.g., soil reflectance (201, 202); LULC (203)) for improving SCORPAN-type models. It is crucial to investigate if RS measures based on spectral analysis can be meaningfully correlated with SOC and other S factors, or to properly represent O factors that are indicative of biological activities or ecosystem types within the modeled pixels. Future work should explore if the combined use of mechanistic functions and RS-derived SCORPAN covariates can reduce the uncertainty associated with simulating soil processes through time (38, 180, 204).

Among the covariates used, precipitation (56% of the cases) and elevation (40%) were most commonly identified as ‘very influential’ ( Figure 4A ) for quantifying agricultural SOC estimation from surface soil layers (< 30 cm) at the regional or broader scale ( Table S2 ). The most commonly identified covariates included precipitation (73%), temperature (62%), elevation (60%), slope (58%), NDVI or other vegetation indices (44%), and land use/cover (38%), topographical wetness index (33%), aspect (23%), multi-resolution valley bottom flatness (19%) when the groups of ‘very influential’ and ‘influential’ were both considered. These influential covariates are mainly comprised of C, R, and O predictors. Key S covariates, including soil type and texture, were less commonly identified as ‘very influential’ or ‘influential’ but had comparable weighted-rankings with other covariates. This is because the S covariates were less frequently chosen for SCORPAN model building ( Table 2 ). Future large spatial scale work should incorporate S covariates more often in order to compare their utilities with other SCORPAN predictors.

According to studies carried out to quantify agricultural SOC from subsurface soil layers (> 30 cm) at the regional or broader scale ( Table S2 ), temperature (48%), elevation (41%), slope (41%), and precipitation (37%) were most frequently identified as ‘very influential’ ( Figure 4B ). By combining ‘very influential’ and ‘influential’ groups together, we identified covariates with good utilities that included slope (78%), precipitation (63%), topographical wetness index (59%), temperature (59%), elevation (59%), aspect (37%), valley depth (37%), and vegetation indices (37%). The relative rankings weighted by the frequency of covariates being used in SCORPAN models emphasized the significance of R covariates involving valley depth, mid-slope position, multi-resolution ridge top flatness, and hill shade. Even though agriculture area is generally located in flat area where the influence of R covariates are likely reduced, it seemed difficult to capture spatial dynamics of subsurface SOC with covariates or measures (e.g. from RS) that focus on surface heterogeneity. Miller et al. (205) pointed out that SOC levels of soil subsurface layers were largely influenced by hydrologic factors rather than O factors, which is consistent with our finding that some of the key C and R covariates that can greatly affect hydrological processes were identified as ‘very influential’, while vegetation indices and land use/cover were more commonly identified as ‘influential’ rather than ‘very influential’ for this category.

For quantifying agricultural SOC estimation from surface soil layers (< 30 cm) at the field or local scale ( Table S2 ), elevation (65%) and vegetation indices (42%) were more frequently identified as the ‘most influential’ covariates than others ( Figure 4C ). The ‘most influential’ together with ‘influential’ groups identified elevation (77%), topographical wetness index (50%), vegetation indices (46%), slope (42%), and multi-resolution valley bottom flatness (31%) with good model utility. The results differed from the ones identified by regional and national scale studies, in that only R and O covariates seem to be most important at this spatial scale. Typical C covariates such as temperature and precipitation were not as frequently identified as influential in these studies. Climate conditions as the first-order control have more predominant influences on SOC at a broader spatial scale, while second-order controls such as microtopographic covariates are better at describing erosion and water flow over small areas and therefore are more influential on SOC distribution at the local scale (206, 207).

Among studies quantifying agricultural SOC estimation from subsurface soil layers (> 30 cm) at the field or local scale ( Table S2 ), elevation (47%) and slope (33%) were most frequently identified as ‘most influential’ ( Figure 4D ), and covariates identified as ‘influential’ or ‘very influential’ most commonly included slope (80%), elevation (73%), topographical wetness index (47%), plan curvature (47%), and profile curvature (40%). It should be noted that the dataset is too small for this category so our results might be biased towards studies with more records than others. However, the observations we reviewed here showed much stronger utility of R than other covariates. This is in line with our regional and national scale results where O covariates showed weaker utility in quantifying agricultural SOC from subsurface than surface soil layers compared to R covariates. It should be noted that the comparison of importance between O and R covariates might differ for regional and national scale studies that involve more heterogeneous ecosystems and vegetation types.

The identification of covariate importance is critical for building effective SCORPAN models at the large spatial extent because use of redundant covariates would cause additional computation costs and introduce sources of uncertainty. However, it is unrealistic to generate a universal set of important covariates because the covariates identified are influenced by the characteristics of the studies (e.g., site characteristics, spatial extent and resolution, soil depth, modeling methods, data sources and interactions among covariates). For example, the R covariates were illustrated to be more influential for SOC estimation compared to S or P factors in higher elevation areas where lower temperature favored the efficiency in SOC stabilization, while the opposite were found in lower elevation regions (208). Viscarra Rossel et al. (209) also reported the significance of regional controls on SOC mapping for large spatial extent and the importance of considering the interactions between C and S, R, and P covariates. These findings call for regionally or locally calibrated SCORPAN models (e.g., separate calibration samples and covariate selection for each predictive point or group of predictions), which, to our knowledge, has not been incorporated to the regional and national scale modeling scheme but should be plausible in the future with the advancement of high-performance computing platforms. Moreover, even though our work points to use of C, R, and O predictors as generally effective SCORPAN covariates for estimating agricultural SOC from the large spatial extent, more efforts are needed to confirm the role of S or P factors regarding whether they contribute to additional predictive power that cannot be achieved with regularly used covariates. It might be the case that S or P factors can be used effectively to delineate representative landscapes (210) for regional or national scale within which the C, R, and O predictors could be applied to build robust SCORPAN models for SOC estimation.