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.

This article proposes an iterative correction method to update soil property maps based on the availability of new soil property measurements. The article is broadly well written and represents a unique and useful contribution to digital soil mapping; however, I have significant reservations about the current state of the methods description and recommend MAJOR revision. I honestly did not read past section 2 after spending > 3 hours trying to understand the methods section because a close reading seemed to reveal inconsistent terminology and only made me more confused. Figure 4 was an attempt to graphically explain the model, but it only confused me more.

Also, I initially understood that this method was only applicable to POLARIS or pHRF predictions, but the author’s response to the comments from another reviewer stated that this method was applicable to ANY soil map. 

I am also deeply concerned about access to the code to reproduce this research. This seems like a method that I would be interested in applying to other soil map products, but the authors only state that the data (not the code?) is available upon request. Why not put this on github?

I'm also confused, does this operate on a cell-by-cell basis or is the distribution adjusted for all grid cells?

I am hopeful that the following comments will be useful to convey areas that I found confusing.

Dear reviewer,

Thank you for your constructive suggestion. We appreciate the time you invested in carefully evaluating the manuscript. Below, we respond to your main concerns and outline how we will revise the manuscript accordingly.

Applicability to Any Soil Map: We will revise the manuscript to state that the proposed iterative residual correction method is not restricted to POLARIS or pHRF products. Instead, it is applicable to other soil map products that provide a spatially continuous (prior) estimate of soil properties within the study area. And we need such information to calculate residuals (difference between new observation and prior estimate of soil property values), since residual is our target variable. The method's applicability to other soil map is inherent in this design. Other soil dataset that provides a continuous baseline estimate can be used as the "prior map" in this process.

Here we emphasize the requirement for using a spatially continuous soil map, because this method does not require prior and new observations to be collocated at the same pixel. In practice, new and legacy soil samples are rarely at identical locations. More likely, at locations where new soil observations are available, there are no collocated legacy data present. The producers of the prior soil map normally use nearby information to interpolate estimate of soil property over this pixel. Hence, the main idea is to update existing soil maps when newly collected data becomes available

Cell-by-Cell vs. Global Adjustment: This method performs cell-by-cell updates and produces region-wide distribution adjustment. The process begins by sampling residuals at grid cells where new soil observations are available. A Random Forest regressor is then trained on these cells to learn the relationship between the residuals and the feature space. Once trained, the regressor interpolates this learned relationship across the entire study region.

As a result, the distribution of soil property values is updated individually at each grid cell through predicted residual corrections. At the same time, because the residual correction is applied consistently across all cells within the study region, these local updates collectively alter the region-wide distribution of soil properties.

Code Sharing: We will update our "Code and Data Availability" statement and already made a public GitHub reposotiry to share the code. Here it the link: https://github.com/emmaxu43/IRC_CA/tree/main

Clarifying Figure 4: We will revise Figure 4 to more clearly distinguish between point-based observations, prior map, components of feature space, predicted residuals, and the interpolated posterior correction. The revised figure 4 will use a more concrete example to illustrate the iterative correction process. Additional details on how the figure will be updated are provided in our point-by-point response to Referee #3. Here, we would like to briefly explain the method:

1. Residual calculation: At locations with new soil observations, residuals are computed as the difference between ground-truth measurements and values from the prior soil map (not restricted to POLARIS or pHRF).
2. Spatial modeling of residuals: A Random Forest regressor is trained using these residuals as the target variable. The feature space includes environmental covariates (e.g., topography, remote sensing data) and soil covariates (e.g., prior soil property values and top-probable values from previous iterations), sampled at the same locations.
3. Map update: The trained model predicts residual corrections for all grid cells in the study area. These predicted residuals are added to the prior map to produce posterior soil property estimates.
4. Iteration: This procedure is repeated until the change in residuals between iterations falls below a user-defined threshold.

Line 14: “Existing soil maps”. Does this mean all soil maps?

Thank you for your suggestion. We agree that the phrase “existing soil maps” is ambiguous. In the revised manuscript, we will clarify this in the abstract by replacing it with “existing probabilistic soil maps — spatially continuous soil property maps that provide a prior estimate of soil properties”. We will provide a definition about this term.

In the context of this study, the prior specifically refers to outputs from the pHRF method. However, the framework is not restricted to pHRF and can be applied to other probabilistic DSM products, provided that a prior prediction surface is available.

Line 27: “prior probabilistic soil property maps”. Does this mean this only applies to POLARIS-like maps? The sentence starting on Line 30. Is this for every pixel in the map?

Thank you for your suggestion. The method is not limited to POLARIS- or pHRF-like products. In this manuscript, pHRF is used as an example to demonstrate the residual correction framework, but the approach itself is more general. Other soil maps can also be used as prior probabilistic soil property maps using this method, please refer to the previous answer.

Regarding line 27: here the “prior probabilistic soil property map” refer to the pHRF-derived soil maps. And the residual correction is applied to every pixel in the map. Residuals are first learned at locations where new soil observations are available, using a Random Forest regressor trained with collocated environmental and soil covariates. The trained model then interpolates residuals across the entire study area, allowing residual corrections to be applied at all pixels within the study region.

Line 73: Latin-hypercube

Thank you for correcting the words. We will make the change in the revised manuscript.

Line 77: Geostatistical models do not heavily rely on expert knowledge. This is an odd sentence.

Thank you for your suggestion. We will replace ‘expert knowledge’ with ‘Traditional geostatistical models often require presumed parameterization and are constrained by stationarity assumptions’ to avoid ambiguity.

We intended to highlight that traditional geostatistical models (such as Kriging) are often constrained by stationarity assumptions and the requirement to select specific semi-variogram models. These models require subjective decisions regarding parameters (such as the nugget, sill, and range); and assume that spatial autocorrelation follows this form in the study area. Without knowledge in the study area, it is difficult to perform these tasks.

Where soil properties exhibit abrupt transitions or complex non-linearities that violate these assumptions, geostatistical models may fail to capture the true spatial heterogeneity. Applying these models without knowledge about the study area can lead to biased estimations.

Line 118. Improves

Thank you for your correction. We will correct it in the revised manuscript.

Line 120. I don't understand this sentence. It does not fit with the rest of the paragraph.

Thank you for your suggestion. In the revised paragraph, we will delete this sentence.

Line 144. Section 2. This section is confusing and I suggest this section is reorganized by first describing the residual correction method, secondly working through a simple example, and thirdly presenting the CA case study and associated data. Step two might not be needed if the authors can more clearly describe the methods, but I would find a clearly worked example helpful for my understanding.

Thank you for your constructive suggestion. This concern is consistent with the comments raised by Referee #3, and we will address both sets of comments through a substantial revision of the methods section, including both text and figures. In the revised manuscript, we will reorganize Section 2 to:

• first describe the residual correction framework in a clear, step-by-step manner;
• include a simple, worked illustrative example to demonstrate how the method operates in practice;
• then present the California case study and associated datasets as a concrete application of the method.

As also outlined in our response to Referee #3, we will update Figure 4 and related descriptions and to visually support the worked example.

Line 228: Figure 3b/Line 233. What does matching prior soils with new profiles actually mean?

Thank you for your suggestion. "Matching" in this context refers to the process of spatially and vertically aligning the new soil profile data with the prior soil map. For each new ground-truth observation (a soil profile at a known location), we extract the corresponding information from the prior map. This involves:

• Spatial Alignment: Using the geographic coordinates (latitude/longitude) of the new profile to locate the exact corresponding pixel in the prior soil map raster.
• Vertical Alignment: Extracting the prior map's predicted value (e.g., soil organic carbon content) for the specific depth interval that matches the depth of the new observation.

This one-to-one pairing allows us to calculate the residual at that specific point in three-dimensional space (latitude, longitude, and depth). We will add relevant explanations in the revised manuscript.

Line 248: classes should be taxonomic classes.

Thank you for your suggestion. We will update the text to specify they are taxonomic classes.

Line 253: soil or environmental covariates? Soil covariates are assumed to be other soil attributes.

Thank you for your suggestion. We will correct this terminology to environmental covariates. These features (e.g., topography, remote sensing indices) are environmental covariates.

Line 254: Please replace the , with and.

Thank you for your suggestion. We will make this change in the revised manuscript.

Line 258: Suggested text: “That prunes less plausible predictions to narrow prediction intervals”

Thank you for your suggestion. We will make this change in the revised manuscript.

Line 273: What does soil data mean?

Thank you for your suggestion. In this context, “soil data” refers to soil property values derived from the prior soil map. Specifically, it denotes the prior estimates at a given pixel and depth interval, including the top-probable soil property values used by the model at that stage of the iteration.

We will revise the sentence in the manuscript to explicitly distinguish between (i) prior soil map values and (ii) observed soil measurements, to avoid confusion.

Line 277: “Additionally, soil covariates are prepared for residual correction at each depth as feature space for predictive models.” This is vague, please explain in more detail. 

Thank you for your suggestion. In the revised manuscript, we will restructure the context to first explain how we construct the feature space before describing its application. We will move the detailed description from Section 2.2.2.1 to an earlier position to provide the definition and components of feature space.

Furthermore, in response to other comments, we will add a new figure to visually illustrate the components of the feature space and how they are updated across the two iteration cycles of our model.

Line 282: Are these randomly selected in the algorithm or just for Fig 4?

Thank you for your suggestion. In the algorithm, the starting layer is randomly selected as part of the iterative process. Regardless of which layer is selected to initiate the process, all layers are eventually updated according to the methodology.

Line 284: what variables make up this feature space?

Thank you for your suggestion. The feature space used for residual prediction is constructed from a combination of variables that capture environmental conditions, vertical soil structure, and prior model predictions. Specifically, it includes:

• Environmental covariates: such as terrain attributes and remote sensing indices;
• Depth information: median values of the soil horizon interval for the layer under modeling;
• Representative soil property values: expected (mean) values of soil properties at each pixel;
• Top-probable soil property values: multiple most likely predictions from the prior iteration;
• Inter-layer residuals: differences between top-probable predicted values across soil layers;
• Weights: associated with top-probable values, representing their probabilities.

These variables are sampled at the locations of the new soil observations and used to train the Random Forest regressor to predict residuals. This construction avoids direct inclusion of the target residuals (residuals) themselves.

And we will revise Section 2 and Figure 4 to better illustrate how these variables are combined and iteratively updated across layers.

Line 284: How are weights updated? I thought these were weights of soil property-soil taxonomic class weights? Does this mean that the soil taxonomic class predictions are also updated?

Thank you for your suggestion. To clarify, the soil taxonomic class predictions are not updated in this method. The “weights” in the feature space refer to the probabilities associated with the top-probable soil property values and are kept fixed throughout the residual correction process. These weights are used to compute representative values but are not themselves updated by the model.

Line 286: Totally lost by this point. What prior values?

Thank you for your suggestion. “Prior values” refer to the soil property estimates used as the baseline at the start of the iterative residual correction. In the first iteration, these correspond to the values predicted by the pHRF method. In subsequent iterations, the “prior values” are updated to reflect the sum of the previous iteration’s residual corrections and the original prior.

We will revise the manuscript and Figure 4 to make this distinction explicit, showing how the baseline values evolve with each iteration, to avoid confusion.

Line 311: What soil observations?

Thank you for your suggestion. "Soil observations" refer to the measured soil property values retrieved from the (new) additional georeferenced soil profiles (WoSIS, SCD, and Salinas Valley measurements) that were not used in creating the prior map.

L314: Feature space is refined… Are the environmental covariates also updated/changed?

Thank you for your suggestion. In response to previous comments, we define the components of feature space include environmental covariates, soil covariates, and weights. The environmental covariates are not updated or changed during the iterative residual correction. Only the soil-related features, such as top-probable values, representative soil property values, and inter-layer residuals, are updated in each iteration. The environmental covariates (e.g., terrain attributes, remote sensing indices) remain fixed and serve as stable predictors for interpolating residuals across the study area. We will revise the manuscript to state that only soil-derived features are iteratively updated, while environmental covariates remain unchanged.

L333: soil property distributions?

Thank you for your suggestion. In this context, “soil property distributions” refers to the cell-by-cell distributions of soil property values across the study region. Convergence is determined based on the median change in soil property values across all pixels: when the median difference between updated residuals and previous residuals falls below a predefined threshold, the iterative updates are considered stabilized. This ensures that further iterations do not largely alter the predictions.

L340: has instead of have.

Thank you for your suggestion. We will correct the word in the revised manuscript.

Fig. 4b. Do these colors mean anything? The rows in 4b seem like they are depth intervals, but these only come from D2? The caption didn’t help much.

Thank you for your suggestion. The colors and rows in Figure 4b are intended to illustrate the multi-layered structure of the feature space and to differentiate among feature columns. The same color indicates the same type of feature, while different rows represent different values of that feature. For example, if we maintain the top-12 probable soil property values at each pixel, the 12 rows correspond to the 12 values for that pixel and depth.

While the example in Figure 4b uses Depth 2 (D2) for illustration purposes, the algorithm updates one depth layer at a time in each iteration, and all soil layers will all undergo this correction process. During an iteration, residual correction is applied only to the “target layer” (for example, D2), while soil covariates from other layers are used to compute inter-layer differences that capture vertical correlations (this feature is illustrated in a previous answer, such as D3-D2 or D1-D2). The algorithm then cycles through all depth layers sequentially, ensuring that each layer is updated. We will update this figure to better illustrate this process.

L343: R3 subscript

Thank you for your suggestion. We will update this correction in the revised manuscript.

L387: So these metrics were calculated as the difference between the values used to iterate... I don't understand this.

Thank you for your suggestion. The metrics were computed using out-of-bag (OOB) samples, which are ground-truth observations not used in training the Random Forest regressor. For each OOB sample, we compare the observed soil property value with the expected value of the updated prediction at that location.

To calculate the updated prediction, at each pixel, we maintain multiple top-probable soil property values from the previous iteration, each with an associated weight. The Random Forest regressor predicts residuals for each of these values in the current iteration. Note the Random Forest regressor is re-trained in each iteration. We then compute the weighted sum of the prior value plus predicted residual to obtain the expected soil property value for that pixel. This expected value is what is compared to the OOB sample to calculate the performance metrics. The reported metrics are computed using all OOB samples and their corresponding expected predictions to summarize the model performance across the study region. We will revise the manuscript to clarify this calculation

L403. Greek rho instead of p?

Thank you for your suggestion. We used the Greek letter rho to represent the Pearson correlation coefficient, rather than the letter P. We will make sure this term is consistently used throughout the revised manuscript.

This paper presents a method to refine predictions of DSM models in a non-parametric manner. I think it may be a significant contribution to DSM and pave the way for promising avenues, such as updating and improving soil maps when newly collected data becomes available. Nonetheless I found the methodology hard to understand (especially how “updated probabilistic maps of soil properties” are obtained and how the “optimization of the feature space” is carried out) and I am concerned it would be hardly reproducible unless the code is made available with a working example. I would also highly recommend the authors to present a minimal example in the paper, on a couple data points.

But I do not think additional experiments or methodological changes are needed: I would therefore recommend to accept the manuscript after some minor revisions.

Comments:

I found the introduction very well written and well supported by extensive references, clearly showing where the current research gaps lie.

Thank you very much for your thoughtful and constructive review, and for recognizing the potential contribution of this work. We agree that the original section 2.2.2 did not sufficiently explain how the posterior soil property maps are obtained. In response, we are pleased to address the points raised:

• Clarification of Methodology: We will revise this section to more clearly explain how the iterative residual correction is implemented and how the updated probabilistic maps are generated. We will also give simple example to explain this method, especially for explaining how we perform iterative update of feature space. We also acknowledge that the term “optimization of feature space” may be misleading, as no explicit optimization objective is used in this process; accordingly, we will replace it with “iterative update of the feature space” throughout the revised manuscript.
• Code Availability: We agree that reproducibility is important. We will share the code via this public GitHub repository: https://github.com/emmaxu43/IRC_CA/tree/main

We will also address these points in more detail in our subsequent, itemized responses to the reviewer comments, where we explain precisely how each suggestion has been addressed. We thank the reviewer for recommending acceptance.

Section 2.1: Your method looks very sensitive to discrepancies and biases across different spatial sources. Therefore I would suggest adding a subsection where you explain more clearly how you ensured that the datasets were consistent with each other (include what you replied to Reviewer #1), and clearly indicate how you divided your data into training/test sets.

Thank you for your suggestion. We will add a subsection as follows to explain how we ensure the consistency of different sources of soil data.  ‘To ensure consistency across the different data sources, we applied several quality control steps. First, we checked the physical plausibility of all soil property values. We defined a valid range with specific minimum and maximum thresholds for each property; any data point falling outside these ranges was considered an error and removed. For soil texture, we required the sum of sand, silt, and clay fractions to equal 100%. If a profile did not meet this compositional constraint, it was excluded.

Furthermore, the datasets are compatible because the WoSIS records for California are largely derived from the NCSS database. Both the SCD and WoSIS datasets follow standardized laboratory protocols, such as those from the Kellogg Soil Survey Laboratory (KSSL). For our own field measurements, we used the Integral Suspension Pressure (ISP+) method to maintain precision for particle size analysis. Vertical soil profiles are also examined. We used equal-area quadratic splines to harmonize all measurements into standardized depth intervals (0-5 cm, 515 cm, 15-30 cm, 30-60 cm, 60-100 cm, and 100-200 cm) for the final analysis.’

Regarding the division of data into training and test sets, we employed an out-of-bag (OOB) sampling. We will explain this in the revised manuscript as well.

Section 2.2.2. In my opinion, this section does not explain how posterior soil properties maps are obtained, but simply how some values are updated based on some available ground-truth data. Until I saw the reply you made to Referee#1, I had not understood either how the method applied to pixels with no ground truth observations. You say that you interpolate the residuals by using a model which learns the correlations between the residuals and other environmental covariates: what model is it? Another pHRF? It should be explained in this section.

Thank you for your suggestions. We will revise the section 2.2.2 to clarify how the posterior maps are obtained. The posterior soil property maps are obtained through the following steps:

• Residual computation: At locations where new soil observations are available, we compute residuals as the difference between the newly observed soil property values and the corresponding values extracted from an existing, spatially continuous prior soil map.
• Residual modeling: A Random Forest (RF) regressor is trained to learn the relationship between these residuals and a set of environmental covariates.
• Residual interpolation and posterior update: The trained RF model is applied over the entire study area to generate a spatially continuous residual correction map. The posterior soil property map is then obtained by adding this residual map to the prior soil property map. Only values of soil property values are updated. Their weights stay the same.

Regarding the applicability of the method to pixels without ground-truth observations, we acknowledge that our earlier response to Referee #1 may have caused confusion, and we will clarify this point in the revised manuscript. In this study, “pixels with no in situ observations” refers to pixels that do not have co-located soil measurements, but there are soil observations existing nearby. A prior estimate of soil property value can be obtained using spatial interpolation for this specific point or pixel. The proposed residual correction method requires the existence of an initial soil property map (the prior) as a model input.

Importantly, our method does not aim to infer soil properties in regions where no soil observations exist anywhere in the study area. Such a scenario falls outside the scope of this work. Instead, the contribution of our approach lies in updating and improving an existing soil property map when new soil observations become available, by learning and performing residual corrections to improve the estimates of soil properties.

We will revise the manuscript to distinguish between (i) pixels without co-located measurements and (ii) regions without any soil observations, and to clarify that the latter case is beyond the scope of this work.

l.273: “By adding these residuals to the prior distributions, the statistical shape of the probability distribution is adjusted (updated property; UP).” Are the residuals really added to prior distributions, or to prior values? From what I understood, residuals are added to prior values.

Thank you for your suggestion. You are correct: residuals are added directly to the prior values at each pixel. At each pixel or sample location, the posterior value is obtained by arithmetic addition of the predicted residual to the prior estimate, while the underlying weights remain unchanged.

We will revise the text to make this distinction explicit. At the same time, we clarify that although residuals are added at the value level, this operation may have a direct consequence on the statistical shape of the resulting probability distribution. Shifting individual values, without modifying the weights, can change the range and spread of the posterior distribution. For example, in Figure 7, the posterior soil maps exhibit an expanded range of soil property values compared to the prior maps: lower values become lower and higher values become higher. This leads to an altered distribution of probability mass across value bins and therefore alters the overall distributional shape. In the revised manuscript, we will update the text to reflect this distinction.

l.277: “Additionally, soil covariates are prepared for residual correction at each depth as feature space for predictive models.” What does this preparation involves? And more generally, I find your use of the term “feature space” quite confusing.

Thank you for your suggestions. We will replace soil covariates to environmental covariates here. We agree that we did not clearly distinguish between static environmental covariates and quantities derived from intermediate model predictions, which led to an ambiguous interpretation of what constitutes the feature space and how it evolves during the iterative procedure.

To address this issue, we will revise Section 2.2.2 by defining the feature space at the beginning of the section. In this work, the feature space is defined as the combination of (i) static environmental covariates inherited from the pHRF method, which remain fixed throughout the procedure, and (ii) prediction-derived features (e.g., current soil property estimates, and inter-layer difference of soil property values; as shows in Figure 4b); and (iii) weights of soil property values that are updated iteratively after each residual correction step.

We will also revise the terminology throughout the section by replacing “optimization of feature space” with “iterative update of feature space”, to more accurately reflect the fact that no explicit optimization objective is defined over the feature space itself.

In addition, we propose to remove the sentence here, as it attempts to mention the preparation of the feature space before Section 2.2.2.1. Presenting it out of order may confuse readers, since the detailed definition and components of the feature space are explained in that subsequent subsection 2.2.2.1.

Section 2.2.2.1 Similarly, I don’t understand what you mean by “optimization of feature space”.

Thank you for your suggestion. By “optimization of feature space,” we were referring to the iterative process of updating and interacting features (static environmental covariates and other updated soil covariates illustrated in the previous answer), as the model progresses.

Rather than using a static set of predictors, the feature space is updated at each iteration to incorporate the most recent updated predictions (UP), allowing the model to capture incremental adjustments and vertical relationships between soil layers that are not reflected in the initial prior map. In addition, we will rephrase this in the manuscript as “iterative update of feature space.” as mentioned in the previous answer.

In the revised manuscript, we will make a plot to compare differences between feature in iteration n and iteration n+1 and show in detail what are the components within the feature space. This figure will be added in the section 2.2.2.1 and provide a clear, illustrated explanation of the feature space and its iterative update, combining both text and figure to improve clarity for the reader.

l.301 I am very confused by these newly introduced terms: what are representative soil property values and how are they obtained? How do they relate to prior model predictions and/or ground-truth data?

Thank you for your suggestions. The representative soil property value at each pixel corresponds to the expected value (mean) of the soil property distribution, representing the current best estimate for that location. For a single pixel, there is only one representative value. In the feature space, for example, if we maintain the top-12 predicted values at a pixel, all 12 entries along this feature column will have the same representative value.

Top-probable values are the multiple most likely predictions generated during the iterative process. Each pixel may have several top-probable values, each associated with a weight reflecting its probability. The representative value can be computed as the weighted sum of these top-probable values.

During each iteration of the Random Forest regression, the model aims to minimize the differences, that is, the difference between the ground-truth observations and the current representative values (predicted residuals plus prior soil property values). Prior predictions are used to compute the representative value at the start of each iteration; for example, in the first iteration, the representative value equals the prior prediction plus the first predicted residual. Subsequent iterations update the representative value by adding newly predicted residuals, progressively refining the soil property estimate at each pixel.

l.303 Aren’t “top-probable soil property values” prior model predictions?

Thank you for your suggestions. While the top-probable soil property values are related to prior predictions, they are not identical. The prior predictions provide the initial estimate of soil properties at each pixel before any residual correction. During the iterative process, the top-probable values are updated predictions that incorporate both the prior and the residuals predicted by the Random Forest regressor at each iteration.

In other words, prior predictions reflect the previous estimate of soil property values; and top-probable soil property values means adjusted estimate of soil property values (residuals plus priors). We will add relevant clarification in the revised manuscript.

l.306 What is the “target layer” ?

Thank you for your suggestions. In our manuscript, the term “target layer” refers to the specific standardized soil depth interval (e.g., 0–5 cm, 5–15 cm) that is currently being processed and updated by the model during a given iteration. In our method, residual corrections are applied depth by depth. During each iteration, residuals are adjusted for only one target layer, while information from other layers is used to provide vertical context.

Because our method processes soil properties depth-by-depth, the "target layer" is the one receiving the residual correction, while "other layers" refer to the other vertical intervals used to provide context of vertical profiles. We will update the text to define this more clearly in the revised manuscript.

What are the dimensions of this “feature space”? And, from my machine-learning experience, the feature space does not include the target variable, only the explanatory features/variables. It seems also very confusing to me that the residuals are part of this feature space.

Thank you for your suggestions. We use D to denote the dimensions of the feature space. The total dimensions is 30. The feature space is composed of the following components:

D_feature = D_{static env covar} + D_depth + D_{mean spv} + D_{top spv} + D_{inter-layer difference} + D_weights

D_{static env covar} = 21

These are the environmental covariates used in the pHRF framework.

D_depth = 1

It is the median values of the soil horizon interval for the soil layer under modeling (current layer).

D_{mean spv} = 1

It is expected value of soil property value (spv) at each pixel.

Dtop spv = 1

It is top-probable soil property value. On each pixel, we can, for example, have top-12 soil property values. And we can flatten all these values and put them within one feature column.

Dinter-layer difference = 5

With six depth intervals, the residuals between the current layer and the other five layers are included to capture vertical variation.

Dweights = 1

The weights associated with each top-probable soil property value.

Summing these components, the total dimension of the feature space is 30.

The feature space does not directly include the target variable for the current layer (i.e., the residuals being predicted). Instead, it includes features derived from prior information and they are not directly collinear with the current target variables, such as:

• Residuals plus prior soil property values (i.e., top-probable soil property values from previous iterations), computed by adding the predicted residual from the previous iteration to the prior prediction.
• Expected value of top-probable soil property values, computed as the weighted sum of the top-probable values at each pixel.
• Differences in top-probable predicted soil property values between layers, calculated as the difference between top-probable values of the current layer and other layers.

These features are not simply linear transformations of the target variable, and they capture incremental adjustments and vertical relationships.

Figure 4 is full of details that I found a bit confusing, e.g. the “prior distribution” at a random pixel, that is unused, uncommented and does not refer to either A or B. The weights, which are a methodological detail that is not properly explained and thus remains unclear to me. I would highly suggest updating this figure using a more concrete example, with the prior and posterior distributions at 2 real points A and B, at two different depths.

Thank you for your constructive suggestions. We agree that Figure 4 in the original manuscript is dense and may be confusing. In the revised manuscript, we will update Figure 4 to make it more concrete and interpretable. Specifically:

• We will show the prior and posterior distributions at two real pixels, A and B, each at two different soil depths, to illustrate how the residual correction updates soil property estimates.
• The weights associated with top-probable values will be explicitly annotated and explained, clarifying their role in computing the representative soil property values.
• We will expand the figure to clearly list the features included in the feature space and illustrate how they are updated across iterations.

Section 2.2.2.3: Using this minimal example, you could show how and at which step this “optimization with constraints” takes place.

Thank you for your suggestion. In the revised manuscript, we will use the minimal example introduced in Section 2.2.2.2 (illustrating residual correction at two pixels and two depths) to explicitly show how and at which step the “update with constraints” is applied. Specifically, we will illustrate how residuals are first added to prior soil property values, how any values exceeding physical bounds are adjusted to the nearest bound, and how remaining residuals are redistributed to maintain consistency.

Table 1 : It would be insightful to know how many ground-truth data points were used to refine predictions for each soil property.

Thank you for your suggestions. We will add this information in the revised manuscript.

Figure 5 is very insightful and well commented. Similarly, the histograms on Figure 7 that clearly show the fact that the distribution is more spread-out than prior predictions.

We appreciate your positive comments. These maps show that the posterior soil property maps, after residual correction, reveal more spatial variability compared to the prior predictions, capturing soil heterogeneity that is smoothed out in the prior maps.

Citation: https://doi.org/10.5194/egusphere-2025-5107-RC3