Abstract: Abstract: Soil plays a very important role in the food chain, and hence is a very important pathway through which humans come into contact with most pollutants. Therefore, there is considerable interest in the best way to monitor the quality of the soil to ensure that it is managed sustainably. However, when the need to monitor the status of soil heavy metals for one area continuously occurs, the sampling and analysis procedures are expensive and time-consuming. Therefore, space-time interpolation is necessary because we can use previous soil sampling points to predict present spatial distribution with fewer soil samples. In this paper, spatio-temporal kriging was utilized to model and predict the spatio-temporal distribution of soil heavy metals. The main objectives of this study were 1) to explore the methods of obtaining an experimental spatio-temporal semivariogram; 2) to fit models for experimental spatio-temporal semivariogram; 3) to perform the algorithm of spatio-temporal kriging interpolation; 4) to evaluate the accuracy and uncertainty of spatio-temporal kriging under the conditions of different neighborhoods; and 5) to predict the spatio-temporal distribution of soil heavy metals of a study area using spatio-temporal kriging. The study area was east of Qingshan district, Wuhan city, Hubei province, China. To monitor the degree of soil contamination, we collected topsoil samples from the study area every year from 2011 to 2014. The number of soil samples from 2011, 2012, 2013, and 2014 were 45, 48, 55, and 48, respectively. The concentrations of cadmium (Cd), copper (Cu), lead (Pb), and zinc (Zn) were analyzed and used as experimental data. For Cd, Cu, and Zn, soil concentrations showed a constant increase from 2010 to 2014. However, the concentrations of Pb showed an increase from 2010 to 2013, followed by a small decrease in 2014. The results of K-S tests showed that Cd, Pb, and Cu did not follow a normal distribution, however, Zn followed a normal distribution. Therefore, the data of Cd, Pb, and Cu were transformed to their common logarithms to achieve a normal distribution. As the results of experimental spatio-temporal semivariance, the T parts of semivariograms for LogCd, LogCu, and LogPb were modeled with a linear model, for Zn was modeled with an exponential model; The S and ST parts of semivariograms for LogCd, LogCu, LogPb, and Zn were modeled with a spherical model. The method of fitting models was the genetic algorithm proposed by the author in 2011. With the results of theoretical variation semivariogram models of LogCd, LogCu, LogPb, and Zn, spatio-temporal kriging were performed. To determine the influence created by the number of neighborhoods, we decided to predict the unmeasured ST point using 4 to 20 of the nearest ST sampling points around the predicted ST site. The results showed that including more neighborhoods could result in less prediction variance. However, more neighborhoods might not produce less RMSE. In addition, ordinary kriging was performed using the same year sampling points while generating a spatial distribution for one year. The results of comparison RMSE generated by ordinary kriging and spatio-temporal kriging showed that spatio-temporal kriging can produce higher prediction accuracy than ordinary kriging. The results of spatio-temporal distribution generally reveal a tendency of Cd, Cu, and Zn concentrations to spread from the south-western part to the whole study area over time, while Pb contamination tends to concentrate mostly on the northern and western parts. The paper showed the computational process of spatio-temporal kriging and its application to soil heavy metals. The results showed that spatio-temporal kriging can improve the prediction accuracy with the help of multi - temporal data.