Abstract: Abstract: The availability of soil organic carbon density (SOCD) information is of great importance for the development of ecological agriculture and the study of global climate change. Compared with traditional laboratory analysis, Visible and near-infrared (VNIR) reflectance spectroscopy has proven to be a rapid, non-destructive and cost effective method for estimating a variety of soil properties. It has got rapid development and has been applied in the field of soil research. In the prediction of soil properties by using spectral information, however, traditional linear regression models often ignore the spatial heterogeneity and dependency of soil, and fail to consider the spatial structure of the error term. With the aims to fill the current gap, a total of 232 topsoil samples were collected in Jianghan Plain with their spectral reflectance and SOCD measured. Partial least squares regression (PLSR), ordinary kriging (OK), co-kriging (COK), and regression kriging (RK) were used to estimate SOCD by using differently pre-treated spectral reflectance. Due to the facts that spectral pretreatments are crucial to reduce the physical variability and particle size effect, and are helpful to remove both additive and multiplicative effects in the spectra, five combinations of spectral pretreatments were utilized while predicting SOCD with PLSR. They were Savitzky-Golay smoothing (S-G), S-G + Multiplicative Scatter Correction (MSC), S-G+Standard Normal Variate (SNV), S-G + first derivative (1st), S-G + second derivative (2nd), S-G+1st+MSC, S-G+2nd+MSC, S-G+1st+SNV, S-G+2nd +SNV. The prediction capabilities of the models were evaluated by R2, root mean squared error (RMSE), and ratio of performance to deviation (RPD). Results showed that the RK approach which utilized soil spectra information outperformed the others, with the highest R2-Pred 0.617 and RPD 1.614, and the lowest RMSEP 0.865 kg/m2. PLSR took the second place with R2-Pred 0.605, RPD 1.523 and RMSEP 0.917 kg/m2, which was also acceptable for SOCD prediction. COK and OK generally failed in the predictions of SOCD, with R2-Pred equaled to 0.007 and 0.004, RPD equaled to 0.903 and 0.874 and RMSEP equaled to 1.547 and 1.597 kg/m2, respectively. Results indicated that the RK model, which considered both the spectral reflectance and the spatial structure of the error term of multivariate linear regression model can improve the prediction accuracy of SOCD. The fundamental reasons could be that soil spectra are comprehensive reflections of soil properties and those environmental factors that influence the formation of soil. Therefore, soil spectra are related with the variation of SOCD, and could be helpful in the prediction of SOCD. Besides, the optimal spectral pretreatment for PLSR modelling of SOCD is the combination of smoothing, first-order derivation and SNV. In summary, soil reflectance spectra in the visible and near-infrared region (350-2 500 nm) could serve as an effective proxy variable for SOCD estimation. Given that soil VNIR reflectance spectra are easy and quick to measure, and the measurement is also environmentally friendly, we would like to argue that soil spectral reflectance could serve as an ideal auxiliary variable for the spatial interpolation of SOCD.