The Blue Nile (BN) region is located in south-eastern Sudan, positioned between longitudes 33° 25' 15.71″ E and 34° 48' 02.18″ E and latitudes 9° 49' 08.52″ N and 12° 59' 43.86″ N (WGS 84 UTM zone 36N). Spanning an area of approximately 38,150.3 km² (Figure 1A), it represents a significant agricultural region within Sudan. Elevation ranges from 370 m to over 1,200 m above sea level (Figure 2A). The predominant landscape comprises flat to gently undulating plains (slopes 0%–5%), with localized steeper terrain (5%–15%) along drainage channels and the Ethiopian escarpment margin. Broad convex landforms are particularly evident in the southern portion of the plain (Figure 2B). Topographical attributes were derived from the 30 m (1 arc-second) Shuttle Radar Topography Mission (SRTM) digital elevation model, which has a reported vertical RMSE of approximately 5-10 m (Golin et al., 2024).

The Blue Nile region features diverse land use/land cover, including horticultural areas, agricultural fields, grasslands, forests, and barren lands (FAO, 2012). The climate is semi-arid tropical (Khir-Eldien and Zahran, 2017), with average annual rainfall ranging from 225 mm to 751 mm (Figure 2C) and maximum annual temperatures between 38.8 °C and 44.5 °C (Figure 2D). Situated within the woodland savannah zone, natural vegetation primarily consists of Acacia spp. tall grass forests, alongside Combretum sp. (Habil), Acacia seyal (Talih), and Anogeissus leiocarpa (Sailack) (FAO, 2012; Elnashi and Ahamed, 2014). Regional soils are predominantly alluvial, heavy cracking clays (Vertisols), with clay content ranging from 40% to 60%, dominated by montmorillonite as the primary clay mineral (Finck, 1961; Fadl, 1971; Ahmad, 1983; Sulieman et al., 2018). According to Soil Taxonomy (Soil Survey Staff, 2014a), these are classified primarily as Typic Haplusterts and Chromic Haplusterts.

The soil dataset was obtained from the Land Evaluation Division of the Agricultural Research Corporation in Wad Medani, Sudan, via the SUSIS portal. This dataset includes 2,235 soil samples collected in previous land-use studies, the spatial distribution of which is illustrated in Figure 1B. The collection consists of 1,980 legacy samples recorded between 1981 and 2020 and 255 contemporary samples collected in 2024. At each sampling location, soil was extracted from up to five depth intervals: 0–30, 30–60, 60–90, 90–120, and 120–150 cm. This process yielded 554, 483, 426, 398, and 374 samples for each respective horizon. SOC was determined using the wet digestion method established by Nelson and Sommers (1996). Observations indicated that SOC content decreased with depth, with mean values of 4.93%, 3.78%, 2.41%, 0.56%, and 0.23% across the five successive layers.

Legacy (1981–2020, n = 1,980) and 2024 field samples (n = 255) were harmonized through (1) Standardization to equivalent soil mass (ESM) to account for bulk density measurement differences (Wendt and Hauser, 2013); (2) Analytical method correction: Walkley-Black SOC values were multiplied by 1.15 to match dry combustion equivalents based on regional validation (n = 45 paired samples); (3) Positional accuracy screening: Legacy points with >100 m GPS uncertainty were excluded; (4) Land-use consistency check: Samples from areas with documented land-use change (urban expansion or irrigation development) were removed; and (5) Depth harmonization: All profiles were standardized to 0–30 cm using equal-area spline functions (Bishop et al., 1999). This reduced the combined dataset from 2,235 to 554 spatially unique, quality-controlled points for modelling.

Available Landsat 9 OLI/TIRS environmental covariates that were used were the Normalized Difference Vegetation Index (NDVI) and Bare Soil Index (BSI). Climatic data included precipitation, temperature, relative humidity, and the Standardized Precipitation–Evapotranspiration Index (SPEI) derived from CHIRPS v2.0 and WorldClim version 2.1, which is a 30-year climatological normal (1994–2024). The Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) was used to calculate the terrain attributes: elevation, slope, aspect, topographic wetness index (TWI), and LS factor. The land use/land cover data, obtained from the FAO SWALIM regional classification (2012), were updated through field verification during the 2024 soil survey.

All covariates for each raster were preprocessed prior to building each model, and all rasters were resampled to a common spatial resolution of 30 m to ensure the variables for spatial consistency(climatic variables, terrain attributes, and spectral indices) were re-sampled by bilinear interpolation, while the categorical LULC layer was re-sampled by the nearest-neighbor method to maintain class integrity. The covariate sources, native spatial resolution, acquisition dates, and resampling processes are detailed in Table 1. Topographic variables were assumed to be static throughout the entire study period; climatic variables were aggregated in 30-year means to describe soil-forming processes; and NDVI data were derived from April 2024 data, as it was available around the same time that soil samples were taken and to reflect antecedent vegetation conditions. Using a multi-year average NDVI instead of the April 2024 NDVI data had a negligible effect on model performance, with the coefficient of determination (R²) varying by less than 0.03, and the models were quite robust to moderate mismatches in temporal and spatial scales. The DMWI modified the so-called catchment measures of topographic wetness by including long-term moisture deficits and the local evaporation demand. It defines a new index, which combines usual catchment factors with a drought correction known as DMWI = TWI × 1 ± SPEI normalised (1)

Where SPEI_normalised is the normalized Standardized Precipitation-Evapotranspiration Index, scaled to represent the degree of drought-driven moisture deficit.

The soil physicochemical analyses were conducted according to international standards, classification using Soil Taxonomy (Soil Survey Staff, 2014a), and some physicochemical analyses, which were estimated based on standard methods. Particularly, the hydrometer method was employed to determine the particle size distribution (Beretta et al., 2014). Soil pH was determined potentiometrically with a pH meter (ORION STAR A211) in soil paste. The electrical conductivity was measured in 1:2.5 suspension extracts using a digital EC meter (YSI model 35) according to standard methods (Sparks et al., 2020). At the same time, SOC was determined using the wet digestion method (Nelson and Sommers, 1996). The variables used to calculate SOC stock were SOC concentration (%), bulk density (Mg m–3), soil thickness (30 cm), and coarse fragments (%).

SOC stock was calculated following the FAO GSOC MRV protocol (FAO, 2018; Yigini et al., 2018) Equation 2: S O C s t o c k s   M g   h a − 1 = S O C D 1000 × D 100 × B D × 100 − C O F R 100 (2)

where SOC conc = organic carbon concentration (%), BD = bulk density of fine earth fraction (g cm⁻³), depth = 30 cm (standardized topsoil depth), CF = coarse fragment volume (%), and 0.1 = unit conversion factor (%→ fraction, g cm⁻³ → Mg m⁻³, cm → m, m² → ha).

SOC stocks (Mg C ha⁻¹) were calculated for each sampling point to a depth of 0–30 cm using the standard equation, where SOC concentration, bulk density (BD), soil depth, and the percentage coarse fragments (>2 mm) were taken into account. The mean bulk density across all sampling points was 1.15 Mg m⁻³, ranging from 1.0 to 1.3 Mg m⁻³, indicating spatial variation in soil compaction and texture in the study area.

For the SOC data modelling, the methodologies established by Mondal et al. (2017), Liu et al. (2022) and Zhu et al. (2022) were employed. The data were initially formatted as a comma-delimited CSV file using Microsoft Excel (Microsoft Corp., Redmond, WA, USA) and converted into a shapefile with QGIS software version 3.44.1, utilizing longitude and latitude coordinates for reference. Only SOC data from the 0–30 cm soil depth layer within the current and legacy soil profiles of the study area were included in the modelling, resulting in a total of 554 SOC data points (Figure 1).

Nine environmental covariates (ECOVs) were selected (Figure 3): precipitation (PR), temperature (T), Normalized Difference Vegetation Index (NDVI), Bare Soil Index (BSI), slope, elevation, aspect, LS-factor, and LULC. NDVI and BSI were calculated from Landsat 9 OLI/TIRS bands Equations 3, 4. The NDVI gives an estimation of vegetation health and range from −1 to +1 (Bangroo et al., 2023). N D V I = NIR − R NIR + R (3) B S I = S W I R 2 + R − N I R + B S W I R 2 + R + N I R + B (4)

where near-infrared (NIR), red (R), blue (B), green (G), SWIR1, and SWIR2 are shortwave infrared bands.

Regression Kriging (RK) was employed to estimate SOC stocks. The dataset (n = 554 samples from the 0–30 cm layer) was partitioned via stratified random sampling into calibration (75%) and validation (25%) subsets. Stepwise Multiple Linear Regression (SMLR) was used to identify significant predictors and remove collinearity. The residuals from the regression were interpolated using ordinary kriging. The final stepwise multiple linear regression (SMLR) model included NDVI, clay content, DMWI, BSI, and precipitation as significant predictors (p < 0.01). Multi-collinearity was not present in the model, with the maximum variance inflation factor (VIF) reaching 3.2, well below the standard threshold of 10. The resulting regression equation was

SOC_stock = 15.4 + 42.1 * NDVI + 0.18 * Clay + 8.7 * DMWI - 15.3 * BSI + 0.012 * Precipitation (adjusted R² = 0.64, F = 98.7, p < 0.001).

The residual variogram was modeled using an exponential function fitted via weighted least squares in the gstat R package (Pebesma et al., 2015). The variogram parameters were a nugget of 12.4 (Mg C ha⁻¹)2, a partial sill of 22.8 (Mg C ha⁻¹)2, and a total sill of 35.2 (Mg C ha⁻¹)2, corresponding to a nugget-to-sill ratio of 0.35 and an effective range of 4.2 km. Analysis indicated no significant anisotropy. A Moran’s I test conducted on the residuals confirmed the absence of significant spatial autocorrelation (I = 0.03, p = 0.42), demonstrating that the model sufficiently accounted for the spatial structure of SOC stocks. The model’s performance was then evaluated. The final prediction combined the deterministic trend and the stochastic residual. Model performance was evaluated using Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and the coefficient of determination (R²).

Model validation employed spatial k-fold cross-validation (k = 5) with a minimum buffer distance of 500 m between training and validation points to ensure spatial independence. Although this buffer (<3 km) is less than the residual variogram range (4.2 km), it provides a conservative minimum separation distance that balances the spatial independence assumption with the spatial availability of soil observations. Although some spatial dependence between the training and validation folds was likely, the spatial CV R2 (0.72) was substantially lower than the random CV R2 (∼0.85), indicating that partitioning the data spatially yields a more conservative and more realistic assessment of accuracy (Brenning, 2012; Roberts et al., 2017). This approach prevents overoptimistic accuracy estimates that arise when nearby samples are split between training and validation sets. The prediction error was calculated quantitatively based on Equation 5. e pred = pred i − obs i (5) where e(pred) is the prediction error, predi is the predicted SOC at validation location i, and obsi is the measured value of SOC at that location (Owusu et al., 2020).

The predictive performance of the regression kriging (RK) model for interpolating SOC was quantitatively evaluated using four statistical metrics, namely (1) root mean square error (RMSE), (2) mean absolute error (MAE), (3) mean square error (MSE), and (4) coefficient of determination (R²), as defined in Equations 6–9, respectively. Notably, the RMSE is particularly valuable as it measures the standard deviation of the prediction residuals, highlighting the extent of discrepancy between the observed and estimated SOC values. R M S E = ∑ y i − y ^ i 2 N − P (6) where ŷ_i is the actual value for the observation, yi is the predicted value, N is the number of observations, and P is the number of parameter estimates, including the constant.

The Mean Absolute Error (MAE) and Mean Squared Error (MSE) are two metrics used to evaluate the accuracy of predictions. MAE calculates the average of the absolute differences between forecasted and actual values, while MSE calculates the average of the squared differences. Both metrics reflect the discrepancies between expected and observed outcomes. A MAE of zero indicates perfect predictions (Robinson and Metternicht, 2006). M A E = ∑ i = 1 n y i − x i n (7) M S E = 1 n ∑ i = 1 n Y i − Y ^ i . 2 (8) where (N) is the total number of observation points, (Yi) is the predicted SOC value, and (y^ I) is the actual or observed SOC value. R 2 = 1 − ∑ i = 1 n p r e d i − o b s ¯ . 2 ∑ i = 1 n o b s i − o b s ¯ . 2 (9) where n is the number of data points at ith location, obsi and predi are observed and predicted SOC values, and o̅̅b̅̅s̅is the mean of observed SOC values. The lower RMSE, MAE, MSE, and higher R2 metrics mean more accurate model performance.

To thoroughly evaluate model performance and address the variability in accuracy estimates stemming from uncertainty, we utilized two approaches: confidence intervals and cross- validation with variance estimation. In our framework, the measure of uncertainty is the kriging standard error, obtained from the kriging variance at the prediction location. This combined method guarantees that the accuracy estimates are both reliable (through CIs) and consistent (through CV) by the established best practices for model evaluation (Bayle et al., 2020; Bates et al., 2024). As an aside, the predictor variable coefficients are not used to propagate the uncertainty in the estimates; only the uncertainty due to the spatial interpolation (kriging) is recorded. Propagating that uncertainty would provide a more comprehensive representation of the total prediction uncertainty and is a possible avenue for future work.

Table 2 shows the descriptive statistics of soil properties and soil environmental variables, the Pearson correlation coefficients among them, and the SOC stock observations from which they were calculated (554 SOC stock observations at the sampling point locations). CEC had the highest positive correlation with SOC (r = 0.91, p < 0.001), indicating the close relationship between organic matter and exchange sites on clay mineral surfaces. Organic matter (OM) showed a perfect correlation with SOC and was calculated directly from SOC using Van Bemmelen’s conversion factor (OM = SOC × 1.724), so OM was not used as an independent variable in the regression modeling to avoid circularity. The clay content had a strong positive correlation (r = 0.84, p < 0.01), as expected, because clay minerals are known to protect organic carbon through physicochemical processes. There was a moderate negative correlation between salinity and EC (r = −0.42, p < 0.05), indicating that the organic matter may be adversely affected by salinity in these Vertisols. The most significant positive associations were with the remote sensing indices, namely NDVI (r = 0.72) and NPP (r = 0.78), while the BSI (r = −0.83) was strongly negative and was effective at discriminating between vegetated and bare soil areas.

Correlation analysis identified key SOC drivers (Table 3) and strong positive correlations (p < 0.01) with NDVI (r = 0.72), clay content (r = 0.84), DMWI (r = 0.78), relative humidity (r = 0.82), and precipitation (r = 0.75). Strong negative correlations occurred with BSI (r = −0.83), sand content (r = −0.75), and temperature (r = −0.72). This negative correlation can primarily be linked to the coarse texture of sandy soils, which results in a lower specific surface area, reduced aggregate stability, and enhanced aeration and drainage (Huang and Hartemink, 2020).

The performance of the Regression Kriging (RK) model was evaluated by comparing measured and predicted SOC stocks across different environmental covariate classes (Table 4). Overall, the model showed good predictive performance, with coefficient of determination (R²) values ranging from 0.78 to 0.97, low bias values (−0.41 to 0.32 Mg C ha⁻¹), and RMSE values between 0.30 and 0.78 Mg C ha⁻¹, indicating generally accurate SOC predictions across the study area. Model performance varied among land use/land cover (LULC) classes. The highest predictive accuracy was achieved in forest areas (R² = 0.95, RMSE = 0.43 Mg C ha⁻¹), followed by cropland (R² = 0.91, RMSE = 0.54 Mg C ha⁻¹). Predictions for pasture and bare land also showed good agreement with measured SOC stocks, with relatively small bias values (−0.10 and −0.14 Mg C ha⁻¹, respectively) and RMSE values of 0.38 and 0.40 Mg C ha⁻¹, suggesting that the RK model effectively captured SOC variability across different land cover types. Prediction accuracy also varied with vegetation density. Areas with NDVI ≥ 0.3 exhibited higher predictive performance (R² = 0.92, RMSE = 0.58 Mg C ha⁻¹) than sparsely vegetated areas (NDVI < 0.3, R² = 0.78, RMSE = 0.62 Mg C ha⁻¹), indicating that vegetation cover contributed to improved model performance. Topographic conditions influenced prediction accuracy. The model performed best on gentle slopes (0°–2°), achieving the highest coefficient of determination (R² = 0.97) with a low RMSE (0.45 Mg C ha⁻¹). Prediction accuracy declined on steeper slopes (≥6°), where R² decreased to 0.82, although the RMSE remained relatively low (0.51 Mg C ha⁻¹). Soil texture strongly influenced model performance. The highest predictive accuracy was observed in areas with clay content ≥ 40% (R² = 0.97, RMSE = 0.47 Mg C ha⁻¹), whereas lower clay content (<25%) produced slightly reduced accuracy (R² = 0.89, RMSE = 0.38 Mg C ha⁻¹). Bias values remained close to zero across all clay content classes, indicating minimal systematic overestimation or underestimation. Climate conditions also affected prediction performance. Areas receiving ≥ 150 mm of annual precipitation exhibited high predictive accuracy (R² = 0.94, RMSE = 0.64 Mg C ha⁻¹), while drier areas receiving < 150 mm had lower predictive performance (R² = 0.79, RMSE = 0.58 Mg C ha⁻¹). Nevertheless, bias values remained small in both precipitation classes (−0.24 to −0.41 Mg C ha⁻¹), demonstrating that the RK model maintained relatively stable predictions across contrasting climatic conditions.

Predicted SOC stocks showed significant spatial heterogeneity across the study area (Figure 4). Stocks were highest in the south and east (35–51 Mg C ha⁻¹), moderate in the central region (20–35 Mg C ha⁻¹), and lowest in the north and west (12–20 Mg C ha⁻¹). This pattern aligns with rainfall gradients and land use distribution. SOC stocks varied markedly across land use types.

The results for the land cover categories are presented in Table 5: Agricultural (irrigated): 51.2 ± 8.4 Mg C ha⁻¹, Forest/woodland: 42.1 ± 9.2 Mg C ha⁻¹, Grassland/pasture: 34.9 ± 7.8 Mg C ha⁻¹, Rain-fed agriculture: 28.5 ± 6.9 Mg C ha⁻¹, and Bare land: 14.2 ± 4.3 Mg C ha⁻¹. The standard deviations and ranges observed in bare land (SD = 4.3) were notably smaller than those in irrigated agriculture and forest/woodland (SD = 8.4–9.2). This suggests a greater homogeneity in SOC stocks within degraded, carbon-limited systems.

Conversely, vegetated and managed land uses displayed more heterogeneous carbon dynamics. These differences are likely attributable to variations in local management practices, soil texture, and microclimatic conditions. The elevated SOC levels in irrigated agriculture are attributed to enhanced biomass productivity and regular inputs of organic matter. The substantial SOC stock in forest and woodland areas aligns with expectations of continuous litter input, reduced soil disturbance, and the accumulation of stable, recalcitrant organic matter characteristic of these ecosystems, as noted by Liao et al. (2024). In contrast, the lower SOC stocks in grassland/pasture and rain-fed agriculture likely stem from moisture limitations and comparatively reduced biomass production. Bare land, lacking vegetative cover, demonstrated the lowest and least variable SOC stock, consistent with minimal carbon inputs and increased susceptibility to erosion and oxidative processes (Zheng et al., 2021).

As shown in Table 6, Regression Kriging (RK) achieved the highest predictive accuracy, outperforming both Multiple Linear Regression (MLR) and Ordinary Kriging (OK), with the highest R² and the lowest RMSE, reducing the RMSE by 26% and 9%, respectively. This demonstrates the advantage of integrating regression modeling with kriging of residuals. The standard TWI in the regression-kriging (RK) model produced cross-validation R2 = 0.72, a realistic estimate of the predictive accuracy at unvisited sites (Supplementary Figures 1A,B). The addition of the NDWI covariate improved the model’s overall predictive accuracy, increasing R2 by 0.08 (from 0.72 to 0.80, when compared to the standard TWI model) (Supplementary Figure 1C). The ordinary Kriging was biased by 0.5, meaning it slightly underestimated the predicted value compared to the actual measurements at the validation points. Spatial cross-validation (5-fold with a 500 m spatial buffer) demonstrated good predictive performance, yielding an R² of 0.72, an RMSE of 8.4 Mg C ha⁻¹ (29.4% of the mean observed SOC stock), a normalized RMSE (nRMSE) of 20.3% (normalized by the observed range), a mean absolute error (MAE) of 6.8 Mg C ha⁻¹, and a bias of −0.8 Mg C ha⁻¹, indicating only a slight underestimation of SOC stocks.