the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 21 Nov 2025
Improvement of Soil Properties Maps using an Iterative Residual Correction Method
Abstract. Accurate mapping of soil properties is vital for many applications including precision irrigation and fertilization. However, existing models for digital soil maps underestimate their spatial variability or prediction uncertainties, which introduces risk for applications including agricultural irrigation and fertilization. This study introduces a hybrid approach that combines prior soil predictions with iterative residual correction to improve soil mapping performance using a Californian case study to demonstrate its application. We first generate prior probabilistic soil property maps using a pruned hierarchical Random Forest (pHRF) method. These prior estimates are then refined by integrating additional soil profile data and iteratively adjusting residuals of distribution of soil properties (differences between observation and prior predictions) pixel by pixel. It gradually adjusts the statistical shape of soil property distributions and incrementally corrects bias of prior knowledge with observed soil information. We evaluated soil mapping over California and at 1-km resolution to test the methodology. For residual correction, we compiled laboratory-measured soil profile data from three primary sources: the World Soil Information Service (WoSIS), the National Soil Characterization Database (SCD), and field measurements conducted by the University of California, Riverside (UCR) and the USDA-ARS United States Salinity Laboratory. From the evaluations, the posterior soil texture predictions show an RMSE of less than 10 %, a 7 % reduction compared to the priors (pHRF-derived soil maps). For soil organic matter (SOM) and oven-dry bulk density (BD), the RMSE also decreased, as the priors initially underestimated their spatial variation. Although posterior SOM and BD predictions were less accurate than other soil properties, this was expected since they are dynamic soil properties and their response to environment and anthropogenic activities is more difficult to simulate. The residual correction also showed reduced uncertainties, as demonstrated by narrower prediction intervals compared to the priors. This method also applied optimization with physical constraints, such as ensuring the bounds of soil property values. This study presents a two-step framework that improves accuracy and reduces uncertainty for DSM applications.
- Preprint
(2849 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (until 19 Jan 2026)
-
RC1: 'Comment on egusphere-2025-5107', Anonymous Referee #1, 07 Dec 2025
reply
-
AC1: 'Reply on RC1', Chengcheng Xu, 22 Dec 2025
reply
This manuscript proposes a method to improve the accuracy of digital soil maps by iteratively applying bias correction across soil depths in a spatial context. While this topic is highly important in DSM studies, several points need further clarification to strengthen the quality of the manuscript.
- The iterations are demonstrated for points with observations; how does this method apply to pixels or points with no prior observations?
Thank you for your suggestion. The iterative residual correction method can be applied to areas with no prior in situ observations. The "prior" refers to initial estimates of soil properties at locations where we have newly collected or gathered soil samples. These priors can be derived from various sources, such as legacy soil maps, global soil property products, or soil maps developed from existing in-situ observations. The residual correction uses additional soil samples to calculate residuals (differences between new observations and prior estimates) at sampled locations. The trained model then learns the relationship between these residuals and soil covariates (such as remote sensing data, topography, and depth information), enabling it to interpolate residuals to pixels without observations. This ‘spatial interpolation’ allows the method to refine predictions of soil properties across the entire study area, not just at sampled locations.
- What are the covariates to determine probabilistic maps of soil class and soil properties?
Thank you for your suggestion. To determine prior probabilistic maps for soil class and properties, we used soil covariates including NOAA GOES-R Land Surface Temperature, the USGS National Land Cover Dataset, ESA Sentinel-1 and -2, and USGS gamma radiation and magnetic anomaly data. More details on these covariates are available in the S3 table of the cited paper (line 246 in the manuscript). The covariates used for the posterior soil property predictions (no soil classes were predicted at this stage) via iterative residual correction are listed as the column names in Figure 4b. And corresponding explanations in text are provided in the section 2.2.2.1.
- Need to clearly mention the division of data for model development and the iterative bias correction. To my understanding now, you used data from Wosis and SCD to develop the prior soil maps and used your own samples for correction. If so, the next question is how do you make sure consistency between Wosis and SCD data when combining those for model development?
Thank you for your suggestion. When develop prior predictions of soil class, we used soil pedons from the National Soil Information System (NASIS) and part of SCD; when develop prior soil properties, we then used a harmonized soil properties database mainly developed using SSURGO. Development of prior soil properties can be found in more detail from this cited paper: Xu, C., Huang, J., Hartemink, A. E., & Chaney, N. W. (2025). Pruned hierarchical Random Forest framework for digital soil mapping: Evaluation using NEON soil properties. Geoderma, 459, 117392. And then, when predicting posterior soil properties, we used SCD (the part that did not used in predicting prior soil class and properties), WoSIS, and soil data that USDA and UCR colleagues collected in California.
It is important to make sure soil data from difference sources consistent with each other. In line 186, we mentioned that the SCD soil data used in the study were collected using standardized laboratory methods. SCD is a subset of the National Cooperative Soil Survey (NCSS) database. And WoSIS within the study area are mainly collected from NCSS database. These soil data were collected using standardized methods, such as the Kellogg Soil Survey Laboratory (KSSL). For each soil property, we also set up plausible ranges of soil property values. If any reported value fell outside its range, we removed it from further procedures. For soil texture, compositional constraints were applied. Additionally, in line 218, we mentioned our use of the Integral Suspension Pressure method (ISP+) for precise particle size analysis of soil and sedimentary materials, following standardized methods, when collecting soil samples in California for our measurements.
- While the models were evaluated using the OOB dataset, I assume that you only develop one model for one soil property at a specific depth. How does this approach guarantee that your model is representative of soil variability across California? Also, how did you divide the data for training and validation?
Thank you for your suggestion. We trained one model per soil property across all soil layers. Our methodology ensures representativeness of soil variability across California through several considerations:
(1) Spatially distributed training and validation data: Our soil samples span different areas within CA (Figure 1). The Random Forest model learns relationships between residuals and soil covariates across different environment, capturing heterogeneous soil patterns across the state.
(2) Basin-based prior predictions: The prior soil property maps were developed using HUC8-basin-level models, which better captured local soil heterogeneity compared to a single statewide model. These basin-level predictions were then merged to create the statewide prior. However, for residual correction, the limited number of additional soil samples (relative to those used for priors) required pooling data across larger areas to achieve sufficient sample size for modeling.
(3) Convergence across iterations: The iterative training process demonstrated that the learned relationship between residuals and soil covariates remained consistent across different iterations, indicating stable and reliable spatial patterns.
(4) Vertical correlation preservation: Training a single model across all soil layers, rather than separate models for each depth, allows the model to learn and preserve vertical correlations within soil profiles.
Regarding data division, validation was performed using soil data that were entirely withheld from the training process, specifically separated based on geographic location. By ensuring that validation data were at distinct geolocations distributed across different parts of California rather than clustered in specific regions, we reduced the risk of data leakage and spatial autocorrelation.
Line 159: Please specify the standardised depth intervals. And briefly explain the method of this harmonisation process.
Thank you for your suggestion. We will cite this reference to our revised manuscript Arrouays, D., McKenzie, N., Hempel, J., de Forges, A. R., & McBratney, A. B. (Eds.). (2014). GlobalSoilMap: basis of the global spatial soil information system. CRC press., at line 159. Following the GlobalSoilMap standards, we harmonized all soil profile data into six fixed depth intervals: 0–5 cm, 5–15 cm, 15–30 cm, 30–60 cm, 60–100 cm, and 100–200 cm. The harmonization was performed using a spline function to interpolate soil property values from the original horizon depths to these standard intervals. It is also mentioned in line 585.
Figure 1: I guess the number of data points is different for each property. It might be more informative if you could mention it in the figure.
Thank you for your suggestion. We will add the following clarification to the Figure 1 caption: "Note that the total number of soil measurements varies by property and decreases with depth, with the surface layers and depths below 1 m generally having fewer observations."
Line 177: what are the interval depths for the Wosis dataset?
Thank you for the suggestion. The WoSIS dataset comprise vertical soil profiles with measurements tied to soil horizons and their ranges differ by location. We harmonized all profiles to the six GlobalSoilMap standard intervals (0–5, 5–15, 15–30, 30–60, 60–100, 100–200 cm) using spline interpolation. We will add relevant explanation similar to line 585 for clarity.
Line 186: It might be good to add a reference for the standardised lab method
Thank you for your suggestion. We will add this reference: Arrouays, D., McKenzie, N., Hempel, J., de Forges, A. R., & McBratney, A. B. (Eds.). (2014). GlobalSoilMap: basis of the global spatial soil information system. CRC press., in the revised manuscript.
Figure 2: In which part of California are these samples located? An inset map would be helpful for this case. Is there any sample from Wosis and SCD located around these areas? If so, what is the closest distance?
Thank you for the suggestion. In the revised manuscript, we will include an inset map showing that field measurements were collected between Salinas and Soledad in the Salinas Valley. To illustrate spatial proximity to WoSIS and SCD sites, we will display our collected samples with distinct marker shapes in Figure 1.
Figure 5: This is difficult to see the difference in performance between soil layers. I suggest plotting the evaluation parameter (R2, RMSE, etc) with the depth as the y-axis for each soil property. Also, why is there a negative value on SOM? What does that mean? Units should be provided on each graph of soil properties.
Thank you for the suggestion. We will include a supplemental figure that plots evaluation metrics (e.g., R2, RMSE) by depth for each soil property to make layer-wise differences clearer.
Regarding SOM, the negative values arise because SOM is plotted on a log scale (as noted in the figure caption); negative values indicate very low SOM content. We will also add units directly to each panel axis in the revised manuscript for clarity.
Line 443: What soil properties is this nRMSE value for?
Thank you for your query. In the revised manuscript, we will change this sentence to ‘The radar plots in Figure 6 illustrate improvements from residual correction for sand, silt, clay, pH, bulk density, and SOM, using three normalized metrics: 1−normalized absolute bias (1−|Bias|), coefficient of determination (R²), and 1−normalized RMSE by ranges of soil variability (1−nRMSE)’ to explain the analyzed soil properties.
Section 3.2: Which dataset that you use to evaluate in this section? Please clarify
Thank you for your suggestion. We used out-of-bag (OOB) samples for evaluation. We also removed OOB samples that share same geolocation with training samples. We will clarify this in the revised manuscript as suggested in a previous comment.
Figure 6: Is this the average value across the soil depths?
Thank you for your query. Those metrics are averaged across all soil depths. We will add this classification to the figure caption and context.
Section 3.3: What does the reduced PIW in this case mean? To my knowledge, reduced PIW means that the prediction interval is narrower, but that does not guarantee for more accurate prediction.
Thank you for your suggestion. We agree that a reduced prediction interval width (PIW) alone does not guarantee improved accuracy, as it primarily reflects uncertainty range. Section 3.3 illustrates that the iterative residual correction largely reduced uncertainty in the posterior prediction of soil properties, as the title of Section 3.3 shows it is a paragraph discussing uncertainty analysis. Together with previous sections, our results can demonstrate that the iterative correction method improves accuracy and lowers uncertainty, as supported by the reduced bias and RMSE listed in previous paragraphs in the manuscript.
Line 583-585: This approach should also be mentioned earlier in the method section
Thank you for the suggestion. We will move this description to Section 2.1.3 for a better flow.
Citation: https://doi.org/10.5194/egusphere-2025-5107-AC1
-
AC1: 'Reply on RC1', Chengcheng Xu, 22 Dec 2025
reply
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 160 | 50 | 22 | 232 | 14 | 14 |
- HTML: 160
- PDF: 50
- XML: 22
- Total: 232
- BibTeX: 14
- EndNote: 14
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
This manuscript proposes a method to improve the accuracy of digital soil maps by iteratively applying bias correction across soil depths in a spatial context. While this topic is highly important in DSM studies, several points need further clarification to strengthen the quality of the manuscript.
Line 159: Please specify the standardised depth intervals. And briefly explain the method of this harmonisation process.
Figure 1: I guess the number of data points is different for each property. It might be more informative if you could mention it in the figure.
Line 177: what are the interval depths for the Wosis dataset?
Line 186: It might be good to add a reference for the standardised lab method
Figure 2: In which part of California are these samples located? An inset map would be helpful for this case. Is there any sample from Wosis and SCD located around these areas? If so, what is the closest distance?
Figure 5: This is difficult to see the difference in performance between soil layers. I suggest plotting the evaluation parameter (R2, RMSE, etc) with the depth as the y-axis for each soil property. Also, why is there a negative value on SOM? What does that mean? Units should be provided on each graph of soil properties.
Line 443: What soil properties is this nRMSE value for?
Section 3.2: Which dataset that you use to evaluate in this section? Please clarify
Figure 6: Is this the average value across the soil depths?
Section 3.3: What does the reduced PIW in this case mean? To my knowledge, reduced PIW means that the prediction interval is narrower, but that does not guarantee for more accurate prediction.
Line 583-585: This approach should also be mentioned earlier in the method section