Bayesian Area-to-Point Regression Kriging for Farm-Scale Soil Carbon Assessment

Most soil management activities are implemented at farm scale, yet digital soil maps are commonly available at regional/ national scale. Disaggregating these regional/national maps to be applicable for farm scale tasks particularly in data poor or limited situations. Although disaggregation is a frequently discussed topic in recent 20 DSM literature, the uncertainty of the disaggregation process is not often discussed. Underestimation of inferential or predictive uncertainty in statistical modelling leads to inaccurate statistical summaries and overconfident decisions. The use of Bayesian inference allows for quantifying the uncertainty associated with the disaggregation process. In this study, a framework of Bayesian Area-to-Point Regression Kriging (ATPRK) is proposed for downscaling soil attributes, in particular, maps of soil organic carbon. Estimation of point support variogram from 25 block supported data, was carried out using the Monte Carlo integration via the Metropolis-Hasting Algorithm. A regional soil carbon map with a resolution of 100 m (block support) was disaggregated to 10 m (point support) information for a farm in the northern NSW, Australia. The derived point support variogram has a higher partial sill and nugget while the range and parameters do not deviate much from the block support data. The disaggregated fine-scale map (point support with a grid spacing of 10m) using Bayesian ATPRK had an 87% concordance 30 correlation with the original coarse scale map. The uncertainty estimates of the disaggregation process were given by a 95% confidence interval (CI) limits. Narrow CI limits indicate, the disaggregation process gives a fair approximation of mean SOC content of the study site. The Bayesian ATPRK approach was compared with dissever; a regression-based disaggregation algorithm. The disaggregated maps generated by dissever had 96% concordance correlation with the coarse-scale map. dissever achieves this higher concordance correlation through 35 an iteration process while Bayesian ATPRK is a one-step process. The two disaggregated products were validated with 127 independent topsoil carbon observations. The validation concordance correlation coefficient for Bayesian ATPRK disaggregation was 23% while downscaled maps generated from dissever had 18% CCC. The advantages and limitations of both disaggregation algorithms are discussed. 40 https://doi.org/10.5194/soil-2019-75 Preprint. Discussion started: 30 January 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
The proliferation of digital soil mapping (DSM) and modelling techniques during the past decade has made available vast amounts of digital soil information (Lagacherie, 2008;Grunwald, 2009;Minasny and McBratney, 2016). For example, the GlobalSoilMap project (Arrouays et al., 2014) aspired towards delivering a global 45 coverage of soil attributes at 100 m resolution for all fivesix standard soil depths. The Soil and Landscape grid of Australia (Grundy et al., 2015), maps of topsoil organic carbon content across Europe (de Brogniez et al., 2015), national digital soil maps of France (Mulder et al., 2016), and carbon map of Nigeria (Akpa et al., 2016) are examples of large extent digital soil mapping (DSM).. However, these national or regional digital maps do not capture the sufficient spatial variability of soil for farm-scale decision making . Yet these 50 maps can be converted into useful farm scale information using a combination of fine-scale environmental covariates and geostatistical methods. This process is called downscaling or disaggregation .
Disaggregation is the process of detailing information collected at a small or coarse scale towards a larger or finer scale (Cheng, 2008). Disaggregating techniques can be classified into to two categories; geostatistical and regression techniques. Area to point kriging (ATPK) (Kyriakidis, 2004) and its variants such as area to point 55 regression kriging (ATPRK) are the most common geostatistical techniques applied to disaggregate soil resource information. Incorporating both point and areal data in the kriging system is such a variant of ATPK presented by Goovaerts (2011). That study demonstrated ATPK to map the variability within soil units and secondly introduced a general formulation for incorporating area-to-area, area-to-point, and point-to-point covariance in the kriging system. Replacing the areal mean with summary statistics in the kriging system is another improved ATPK 60 technique suggested by Orton et al. (2014). Román Dobarco et al. (2016) extended this approach using summary statistics in the ATPRK system. Further, a Bayesian regression kriging approachmodel was presented by Ramos et al. (2013). The model used aThey introduced Gaussian process model, and it is as a novel image fusion technique for panchromatic sharpening of low-resolution satellite images using. They used images from unmanned aerial vehicle flown over agricultural fields to demonstrate their model. 65 Regression approaches are also a popular approach to downscale soil attributes. For example, Malone et al. (2012), formalised a statistical downscaling algorithm; dissever based on the Liu and Pu (2008) algorithm for downscaling continuous earth resources with a case study of downscaling soil organic carbon content. This method was refined by Roudier et al. (2017), allowing the user to choose between a range of regression models apart from the originally contained Generalised additive models (GAM). "DSMART" (Disaggregation and 70 Harmonisations of Soil Map Units through Resampled Classification Trees) is another popular regression-based approach, proposed by Odgers et al. (2014) to downscale soil map units to soil series information. However it is different to the discussed method since there is no fine gridding involved.
Disaggregation is a change of support problem (Cressie, 1991;Gotway and Young, 2002). Spatial support is the geometrical size, shape orientation spatial unit of an observation or a prediction. In the 75 disaggregation process, when the support changes the associated statistical and spatial properties of the variable also change, adding a major source of uncertainty (Truong et al., 2014). ATPK is a theoretically-based method that uses the classical kriging principles to predict at point support using block support information (Kyriakidis, 2004). Incorporating fine scale covariate data in the kriging system (ATPRK) is superior to ATPK as it allows for within class variability, spatial autocorrelation and also ancillary covariates representing "scorpan" factors 80 (McBratney et al., 2003) represented by ancillary covariates ; (Kerry et al., 2012).
In disaggregating soil information, quantifying the uncertainties of this process is not often discussed.
Underestimation of inferential or predictive uncertainty in statistical modelling can lead to inaccurate statistical summaries and overconfident decisions (Draper, 1995). Therefore, it is important to account for the underlingunderlying parameter uncertainty associated with the disaggregation process. The use of Bayesian 85 inference in ATPRK system enables quantifying the uncertainty associated with inferring point support parameters from the block supported system, in addition to estimating the product uncertainty (Diggle et al., 1998).
To the best of our knowledge, Bayesian ATPRK has not been used for downscaling soil information. In this paper, we present a Bayesian ATPRK approach for downscaling soil attributes. Specifically, the study focuses on downscaling regional soil organic carbon (SOC) maps to high-resolution farm scale maps. Soil organic 90 carbonSOC is arguably a key soil property due to its conferring benefits to the soil physical and chemical properties and potential to store atmospheric carbon. Accurate estimations of SOC stocks at farm scale are important for precision agricultural needs. In addition, the fine-scale maps can be used for designing sampling strategies for quantification of SOC stocks (de Gruijter et al., 2016) required by carbon sequestration auditing process such as carbon farming initiatives in Australia . 95 A regional soil carbon map with a resolution of 100 m (block support) was disaggregated into 10 m (point support) information for a farm in the Spring Ridge region in the northern NSW. Application of the Bayesian ATPRK approach is presented and demonstrated within the study. In addition, a comparative analysis with the dissever algorithm, a regression-based downscaling method, is performed to assess the performance of Bayesian ATPRK in the context of spatial downscaling of digital soil carbon mapping. 100

Theoretical background
Assuming we are observing a continuous Gaussian process, we denote the process as S(p) for locations p ∈ D, where D is the region of interest. For block support (B) data, where B ⊂ D, we presume that the observations are as block averages: 105 where |B| denotes the area of B.
The underlying Gaussian process S can be expressed as where is the design matrix of the trend function and β vector of coefficients, and ε is the stochastic residuals with mean of zero and covariance function C ( ), where represent the parameters of the covariance function.
Since the block supported observations; and point support observations; are jointly Gaussian, the prediction problem can be introduced as: 115 where is the variance-covariance matrix of , the variance-covariance matrix of , and are the variance-covariance matrices between and and vice versa. is the design matrix of for fine-scale covariates of . 120 Since the joint distribution is Gaussian the optimal predictions at unsampled point locations can be given as, where is the vector of predicted soil carbon at N pointpoints. 125 And the variance-covariance matrix for the prediction error is, This shows that the inference of from is straightforward and resembles a simple kriging formulation (Cressie and Wikle, 2011). However, the estimation of point support variogram from the block support data is not a trivial 130 task. Usually, the estimation is done via deregulation or deconvolution of the empirical variogram. Goovaerts

Bayesian inference
The likelihood function of Gaussian process S can be given as ( )|β, ) = N β, ( ) . is a function of the Euclidean distance h, with the vector of parameters . Then the likelihood of predictive distribution at point locations can be given as, Where | , β, is distributed as Each entry of the distribution can be estimated through a Monte Carlo integration. Then we can replace Eq. 3 with 150 ( , | β , ), where 'hat' denotes a Monte Carlo integration. Then it is apparent that, To predict at we require f , , which is given by Eq. 3. Using ( , | β , ), and using Eq. 8 we can obtain , | β, to sample . 155 Then given the priors [β, ], the Bayesian model is specified, and a Monte Carlo fitting procedure can be carried out to maximise following log likelihood function. A complete description of deducing the Bayesian area to point inference can be found in Gelfand et al. (2001).
Where is block to block covariance and can be approximately calculated as 165 Where i , j index are block supported observations, k, l index discretised points within a block and K is the number of discretisation points per block (Truong et al., 2014). c , c are the partial sill and the nugget component of the semivariogram; γ .

Study area and coarse scale soil carbon map
We illustrate the methodology by disaggregating a regional topsoil organic carbon map of the University of Sydney E. J. Holtsbaum Agricultural Research Station (Fig. 1). The landholding, also known as "Nowley", is located in the Spring Ridge district of the central/north west slopes region of NSW, Australia. Its area is 175 approximately 2300ha and is managed as a mixed farming enterprise of cattle grazing throughout the year, and wheat, barley, and canola in winter, and sorghum and sunflower in the summer. Cropping is performed on soils derived from basalt parent materials. The average annual rainfall is 600mm.
The available map is a New South Wales (NSW) topsoil carbon map, which was generated using a local regression kriging technique (Somarathna et al., 2016). These maps are analogues to the nationally available soil 180 carbon maps of Australia which are available in separate layers corresponding to the depth intervals definition in the GlobalSoilMap product specifications (Arrouays et al., 2014). The standard soil depth layers are 0-5cm, 5-15cm, 15-30cm, 30-60cm, and 60-100cm.
This NSW soil carbon map has grid spacing 100m x 100m, with block support of the same size. Hence, map spatial resolution and the support are equal, which in turn mean that the pixel values that make up the map 185 represent a predicted mean for the area of each pixel. Commonly the pixel values of a digital soil map are on point support, with the pixel value representing the value at a point usually at the centre of each pixel (Malone et al., 2013). Therefore, block kriging was used convert the maps to point supported maps before the analysis.
The topsoil of the study area was represented by the depth interval of 0-7.5 cm, as this was the depth interval used in previous work conducted at this site for soil carbon auditing (de Gruijter et al., 2016). The 190 respective coarse (NSW) map for the 0-7.5cm depth interval was predicted by using mass-preserving splines (Malone et al., 2009) of the 0-5cm and 5-15cm standardised maps. This map will be hereafter referred to as the coarse scale map. It is this map that subsequent descriptions of spatial downscaling will refer to. The experimental variogram of the block supported coarse scale map is given in Fig. 2.

Environmental covariates
Elevation, topographic wetness index (TWI), Gamma radiometric Thorium, gamma radiometric potassium, Landsat band 4 and Landsat band 7 were used as environmental covariates in the spatial downscaling model.
Landsat data has a resolution of 30m. Since the coarse map was disaggregated into 10m resolution, the Landsat data were pan-sharpened using the panchromatic band (15m), and then resampled to 10m. 200 The elevation data had been collected during a ground survey using an all-terrain vehicle at with 20m swath. Then the data were processed using a local block kriging approach to produce an elevation data at a 10m resolution. Based on the developed digital elevation model (DEM,), the Topographic Wetness IndexTWI was derived. The gamma radiometric data was derived from aerial survey (Minty et al., 2009) which were resampled to 10m. A more detailed description of data collection and post-processing can be found in Malone et al. (2017). 205

Model parameter inferencing using the MCMC simulation
Suppose our observed soil carbon data (S) follows a second-order stationary Gaussian process which is 210 characterised by the mean (trend ) and the variogram where ~ N( , γ ) where γ = (h, ) ; is a function of Euclidian distance h and vector of parameters of Matérn variogram (Eq. 14). The block averages have the same spatial mean as but a different spatial structure (Gotway and Young, 2002). The block support covariance is also a function of h with parameter vector of but has a different support.
where C is the covariance between observation i and j, h represents the separation distance between i and j, δ denotes the Kronecker delta (δ =1 if i=j and δ =0 when I ≠ j), c + c signifies the sill variance, Κ is the modified Bessel function of the second kind of order υ. Γ is the gamma function, r denotes the distance or 'range' parameter and υ is the spatial 'smoothness'. Calculating requires discretisation and it was done using a regular 220 grid pattern, with equal 10m spacing between discretising points within a block. Then was calculated according to Eq. 11.
The priors were given as informative priors and considered as normally distributed. The coarse scale soil carbon map was regressed against the selected covariates to estimate the approximate prior values for the coefficient of the linear trend (β). Gaussian priors were assigned for both spatial trend (β) and covariance component ( ) of the spatial model. MCMC integration was carried out using the Metropolis Hasting algorithm.
Metropolis Hasting is a computationally efficient algorithm that allows sampling from conditional distributions (Diggle et al., 1998). Calibration of the model was done using two parallel chains of MCMC with 100 000 iterations. 230

Spatial prediction
The mean, 95% upper and lower CI limits of spatial model parameters were calculated using the values of the last 10 000 MCMC iterations. These calibrated parameters were used to predict SOC content at the fine scale. The predictions were made on to the 10m x 10m points based on Eq. 4. 235

Disaggregation using dissever
The proposed Bayesian model also compared with the dissever regression -based disaggregation approach.
dissever is a regression-based downscaling approached proposed by Malone et al. (2012). dissever operates in two steps; initialisation and iteration. At initialisation step, the coarse map is resampled using nearest neighbour 240 resampling technique to match the resolution of available fine scale covariates. Then the resampled map is regressed against the fine scale covariates. Next, the regressed fine map is upscaled through averaging, and then the mass balance deviation (deviation factor; DF) from the coarse map is calculated. The regressed fine map is corrected using the DF, and then the algorithm progresses to the iteration step. In the iteration step, regression and upscaling continue until the mass balance deviation between the upscaled map and the original map is less than a 245 set minimum (threshold). The following flow diagram (Fig. 3) is an illustrative description of the dissever algorithm. This differs from Bayesian ATPRK since it is not focused to address the change of support problem and it is a fine gridding exercise. A detailed and comprehensive account of dissever can be found in Malone et al. (2012) and Malone et al. (2017).

250
Following the procedure mentioned above, the coarse scale carbon map of the area was disaggregated 255 into 10m using dissever. The environmental covariates used in the Bayesian approach were similarly made available for dissever.
Where ρ is the correlation coefficient, σ , σ are the variances of observed (x) and predicted (y) values and μ , μ are the respective means. CCC is scaled between -1 and 1, the latter implying perfect agreement and the former implying perfect reverse agreement. 280 MSE is given by where is the observed values of SOC at independent sampling points, predicted SOC valevalues at the fine scale and n′ number of independent SOC samples at a fine scale. MSE is an indication of the accuracy of 285 predictions, while CCC indicates both accuracy and precision of the disaggregation exercise.

Results and Discussion
The use of LMMlinear mixed model (LMM) as the spatial prediction model allows integration of fine-scale 290 covariate information and spatial correlation in a single model. Unlike regression kriging in which regression and residual kriging are used in two phases of modelling, LMM is a one-step process. We used Bayesian inference techniques to estimate parameters of LMM to deliver uncertainty estimates of model parameters. This section provides a description of posterior parameters of the spatial model followed by a comparison of spatial disaggregated maps using Bayesian ATPRK and Disseverdissever. 295  values of the range parameter are more-or-less similar to the prior values. However, the partial sill and nugget parameters are different from their priors. The posterior partial sill is higher than the prior sill, which agrees with the fact that the coarse scale variability of soil carbon is lower compared to the finer scale (Fig. 5b).  the Bayesian technique provides a range/interval of the most probable value with a 95% confidence. These 320 estimates are more reliable as they are generated using numerous repeated samplings. Accuracy and confidence in estimates are important parameters for planning, designing, and implementation of soil management. Therefore, Bayesian ATPRK provides a comprehensive output of soil carbon estimates at the farm scale.

Comparison of Bayesian ATPRK and Dissever disaggregation methods
The disaggregated maps from both techniques were upscaled using a spatial overlay to assess the accuracy of disaggregation. CCC values were calculated for both techniques. Fig. 7 shows the agreement between the cell 330 values of coarse map and the upscaled disaggregated maps.
Concordance correlation between original coarse values and upscaled product were 87% and 96% for 335 Bayesian ATPRK and dissever, respectively. dissever iterates until the agreement between upscaled and original coarse cell reaches a set minimum, while Bayesian ATPRK output is a result of one simulation considering both the correlation with the covariates and covariance structure of the residuals. This explains the stronger concordance associated with dissever. Fig.8 illustrates the comparison between dissever and Bayesian ATPRK disaggregated product with the 340 coarse scale source map. A closer examination reveals that Bayesian ATPRK product resembles the spatial pattern of SOC more closely than the dissever product while dissever seemingly over predicts the SOC content.
Conversely, dissever is computationally more efficient compared to Bayesian ATPRK. Calculating the block to point and block to block distance matrixes is a computational expensive procedure. The whole study area cannot be computed in a single run using a desktop computer. Therefore, the study area has to be split into tiles, and 345 predictions were carried out for each of the tiles. This resulted in abrupt changes between two tiles in the Bayesian ATPRK disaggregated map. The use of parallel computing, high performing computing can improve the processing time.

Validation of Disaggregated maps
Independent validation results between disaggregated maps using dissever and Bayesian ATPRK indicate that Bayesian ATPRK technique is slightly more accurate compared to dissever produced output in general. 355 CCC value of Bayesian ATPRK was 23% while that of dissever was 18%. Greater CCC indicates more accurate and precise predictions. Mean squared error (MSE) %) of Bayesian ATPRK product was 0.45 while MSE for dissever was 0.49, which indicates that Bayesian ATPRK disaggregation is slightly more accurate than dissever disaggregation.
Although the validation CCC is low, it still is a sufficient approximation in terms of SOC as the SOC 360 prediction accuracy generally remains low. The coarse map is a product of regression for the whole region, which also has a lower concordance correlation with the sampled data. Through the disaggregation process, the errors tend to propagate widening the gap between the actual and predicted SOC content.
Overall, we have demonstrated the use Bayesian ATPRK for downscaling coarse scale rasters to more informative fine scale rasters using available fine scale covariates. The study proves that Bayesian ATPRK 365 provides more accurate results with uncertainty estimates of the disaggregation process. There are few limitations Formatted: Indent: First line: 1.27 cm of this study. Firstly, we did not deal with the uncertainty from all source of measurement errors. The coarse scale map always caries a certain degree of error/uncertainty. Also, there is an uncertainty associated with the splining of the maps. Both of these uncertainties can be treated as measurement errors and filtered by adding as a noise variance along the diagonal of the covariance matrix . Another way of correcting the bias is to use the 370 measured data in the spatial model to reflect the actual ground conditions. Malone et al. 2017 demonstrated that the inclusion of the ground observations improves the accuracy of the disaggregated maps. A further limitation of our method is that the disaggregated map produced from Bayesian ATPRK has a block effect due to the final map is a collection of tiles produced separately due to high computing demands of the model. This titling effect can be eliminated by using high performing and parallel computing solutions to generate the fine scale image at 375 one go. Also, low rank-representations (Wikle, 2010) of the covariance can also be used as possible solutions to computational burden of the process.
Furthermore, although we have used the proposed Bayesian ATPRK in DSM domain, this method can also be used to downscale coarse scale satellite data. Even though there is a wide range of fine resolution satellite images available in present context, there is still sparsity of fine resolution data for certain parts of the world. 380 Spatial downscaling is a solution to this data sparsity. There is a range of disaggregation/ pan sharpening techniques available in the literature. However, very few techniques can tackle the uncertainty of source images and the uncertainty of the disaggregation process. Unlike other methods, Bayesian ATPRK can explicitly help to tackle observation errors of the satellite images.

Conclusions and Recommendations
The study demonstrates that Bayesian techniques can be effectively used to address the change of support problem.
Monte Carlo integration provides the estimates of point supported covariance structure using block supported data along with uncertainty estimates of parameter inference. This study used informative priors of covariance 390 parameters as the soil carbon data were available for the study site. In the absence of prior information, priors can be given as non-informative priors.
Bayesian ATPRK provides an effective way of disaggregating coarse scale SOC maps. The output of this study facilitates policy making and soil management activities by providing the maps with confidence limits which lead to wider understanding of the SOC content of the area. However, the method is computationally 395 expensiveintensive, and the use of high performing computing techniques is recommended for faster and smoother predictions.    We would like to thank the reviewers for providing constructive comments and efforts towards improving this manuscript. We highlight the reviewers' comments in red and our responses are given in black.

Reviewer 1 comments and Responses
General comments: This manuscript proposed the use of Bayesian Area-to-Point Regression Kriging (ATPRK) for spatial disaggregation of regional scale digital soil map to farm scale. This method is able to provide the 705 uncertainty estimate associated with the disaggregation process. Throughout it had lower concordance correlation with the coarse map than dissever algorithm, the independent data showed that Bayesian ATPRK had a higher concordance correlation than dissever. Generally, the proposed method is interesting as it can provide the uncertainty related to the disaggregation process. However, from my point of view, the major uncertainty of disaggregated map comes from the original map, thus it is more important to incorporate the uncertainty of the 710 original map into disaggregation. Is the proposed method able to integrate this? Another limitation is the Results and Discussion section, where authors focused mostly on the results and rarely provided a more general discussion. I look forward to seeing the feedback from the authors on these two issues We have identified the uncertainty of the original should also be considered in the disaggregation process. We 715 have mentioned this in our conclusions section as a recommendation. Our main focus of this paper was to estimate the uncertainty of the disaggregation process which is overly neglected in the current literature. However, the proposed Bayesian method is capable of handling the uncertainty of the original map by incorporating the uncertainty as a noise variance in the covariance matrix. CBB. This enables filtering the noise. Also, we have added more detailed discussion about filtering measurement errors. (lines 341 -348) 720 Another limitation is the Results and Discussion section, where authors focused mostly on the results and rarely provided a more general discussion.
We agree with the reviewer that the discussion section should be little more generalised. Therefore, a more generalised discussion was added at the of the results and discussion section (lines 341 -361).

725
Specific comments: Line 21: Why only mentioned underestimation? How about overestimation?
Here we are referring to the fact that in DSM literature the uncertainty estimates of the disaggregation process are often neglected.
Line 35: What is the difference in computing time? I know Bayesian ATPRK is a one-step process, and it saves computing time? Because sometimes an iteration process can also be efficient. 730 The difference in computing time is approximately 3 hrs Line 45: It is not "five standard soil depths" but "six standard depth intervals".
The typo was corrected Lines 71-72: I do not agree with the statement here, DSMART can also produce maps in a fine grid.
Here, we are referring to the fine riding process where we disaggregated a coarse grid cell (100m) into many (10 735 x 10) fine grid cells. While the rasters created by DSMART based on the most probable class assigned after generating realisations from classification trees. The classification trees are built using the values taken from random samples from each soil class polygon from the legacy maps.
Line 123: What do you mean "at N point"? Line 136: Since you listed many methods for deconvolution of the empirical variogram, it is better to show here why you chose Bayesian estimation. 740 We have included the objective of choosing Bayesian techniques to infer the point support variogram. (line 137-

138)
Lines 187-189: Do you think the proposed framework can integrate the uncertainty from the fitted masspreserving splines as well as from the uncertainty estimates associated to these two SOC maps (0-5, 5-15 cm)?
Yes, the uncertainty from all sources can be accounted using the proposed Bayesian techniques. This is also 745 included in the added new section to the results and discussion section.
I understand that the Bayesian ATPRK can provide the uncertainty in the step of disaggregation, however, I think the major uncertainty comes from the original map, which should be not ignored in disaggregation.
We have identified this as a limitation of the study, and it was further elaborated in this revised version (line 341- There were two reasons behind choosing 10m as the fine resolution. One is the resolution of available data at fine scale and secondly, we believe that 10m maps are informative enough for supporting decision making process at 755 fine scale. In anotherword,10m maps would capture most of the spatial variability of SOC and therefore we selected 10m as our fine scale's resolution.
Line 268: I suggest adding PICP to assess the uncertainty.
We present CI estimates of model parameters and maps. CI values provide an insight to the uncertainty estimates. Therefore, we decided not to include PICPs 760 Lines 341-342: reduce the processing time? However, parallel computing can still not solve abrupt changes between two tiles. Maybe you need a laptop with larger RAM?
Yes, it needs a bigger RAM to disaggregate the whole study area at once.
Line 579: Figure 6 still showed an obvious block effect. In order to compare the original SOC map (Figure 1), please use the same legend and then we can tell the difference between them. 765 Figure 1 was changed to the same legend.
Line 595: The map from Disserver looks smoother.
Disserver disaggregate the whole study area at once resulting a smoother map. Therefore, the maps are smoother.
Technical corrections Line 21: DSM has not been defined in the previous texts. Line 48: DSM has been defined before. 770 Respective corrections were made.
Line 91: SOC has not been defined yet.

SOC defined
Line 200: C2 SOILD Interactive comment Printer-friendly version Discussion paper DEM has not been defined yet and Topographic Wetness Index can be replaced by TWI here. 775 DEM was defined and the Topographic Wetness Index was replaced by TWI Line 209: and should not be in italic. Please have a careful check of all the parameters which should be in italic.

Respective corrections were made
Line 280: predicted SOC values. 780 Typo was corrected Line 285: LMM has not been defined.

LMM was defined
Line 290: In other texts, dissever is used. Please make it consistent.

Reviewer 2 comments and responses
The manuscript is well-written and follows up from earlier work of Malone et al. 2017 on disaggregation of coarse scale SOC maps. This is particular important for facilitating management at the field/farm scale. Although I am not a geostatistical expert, the data anlysis is sound and the results are convincing, although not much different from the ones obtained by the 'dissever' approach. 795 Lines 19-20 Is not there a verb missing in the sentence starting with 'Disaggregating' ?
The typo was corrected Line 21 Define 'DSM' at first time use.

DSM was defined
Line 52 I appreciate that you try to use the original meaning of scale 'large' and 'small', but the word 'or' might 800 create confusion. How about replacing 'or' by 'i.e.'? Line 53 ..into two categories. . .

Typo was corrected
Line 62 Please rephrase this sentence and break it up. e.g. They demonstrated that a Gaussian process model can be used as a novel image fusion technique. They used images from unmanned aerial vehicle flown over 805 agricultural fields for panchromatic sharpening of low-resolution satellite images.
The sentence was rephrased.
All the following typos were corrected.
Line 65 Avoid using 'approach' twice in one sentence.