Three-dimensional geostatistical modeling of soil organic carbon: a case study in the Qilian Mountains, China.

For mapping soil properties in three dimensions the simplest option is to choose a series of depth intervals, and to calibrate a two-dimensional (2-D) model for each interval. The alternative is to calibrate a full three dimensional (3-D) model that describes the variation in lateral and vertical direction. In 3-D modeling we must anticipate possible changes with depth of the effects of environmental covariates on the soil property of interest. This can be achieved by including interactions between the environmental covariates and depth. Also we must anticipate possible non-stationarity of the residual variance with depth. This can be achieved by fitting a 3-D correlation function, and multiplying the correlation between two points by the residual standard deviations at these two points that are a function of depth. In this paper various 3-D models of the natural logarithms of SOC are compared with 2-D depth-interval specific models. Five environmental covariates are used as predictors in modeling the lateral trend. In the 3-D models also depth was used as a predictor, either categorical, with categories equal to the depth intervals (3-Dcat), or continuous (3-Dcon). The covariance of the residuals in 3-D is modeled by a sum-metric covariance function. Both stationary and non-stationary variance models are fitted. In the non-stationary variance models the residual standard deviations are modeled either as a stepwise function or as a linear function of depth. In the 2-D models the regression coefficients differed largely between the depth intervals. In the 3-Dcat model extreme values for the regression coefficients were leveled out, and in the 3-Dcon model only the coefficients of NDVI and aspect changed with depth. The 3-Dcon model with a residual standard deviation that is a stepwise function of depth had the largest residual log-likelihood and smallest AIC among all 3-D models. Based on the cross-validation root mean squared error (RMSE) there was no single best model. Based on the mean and median of the standardized squared error (MSSE, MedSSE) the 2-D models outperformed all 3-D models. Overestimation of the prediction error variance by the kriging variance was less strong with the non-stationary variance models compared to the stationary variance models. 3-D modeling is required for realistic geostatistical simulation in spatial uncertainty analyses.