Abstract Exposure to Lead (Pb) from soil and dust in urban areas continues to threaten human health especially that of children, who are more likely to ingest a higher Pb fraction from contaminated soil. Geostatistical approaches were evaluated for the prediction of soil Pb concentrations. Ordinary kriging (OK) and co-kriging (CoK) were used to predict spatial distributions of both soil Pb concentrations and bioaccessibility. The relative accuracies of the developed predictive models were assessed using of coefficients of determination ( R 2 ) and the root mean square error (RMSE) based on the cross validation method. The accuracy of the final model was evaluated via comparison between predicted and traditional soil measurements. While OK models ( n = 73 ) were unreliable for predicting a Pb concentration surface, correlation between Pb and Zn was a useful mechanism for obtaining better predictions of soil Pb using CoK. The CoK model of log transformed Pb with Zn resulted in the best fitted model ( R 2 = 0.63 ). The percentage relative improvement (RI) for this model was 39% which suggested relatively reliable prediction accuracy. A probability kriging (PK) surface was used to describe the probability of bioaccessibility exceeding a threshold of 17% as an indicator for potential human risk. The areas with highest probability of exceeding the threshold were in agreement with previous risk area divisions related to blood Pb (BPb) levels for children under 5 years of age. The results of this research confirmed that geostatistical methods had the ability to rapidly estimate soil Pb distributions for environmental health risk assessment.