Evaluating of geostatistical approaches in monitoring of spatial variability of some chemical contaminants such as agricultural phosphorus (P) will provide valuable data for large agricultural areas. In this study, performance of varied GIS based geostatistical interpolation methods were tested in assessing of site specific P variability on the apple orchard. For this aim, soil samples were systematically collected from the agricultural apple area using the grid sampling system. The samples were taken at two depths (025 cm and 25-50 cm), the distance on the Y direction is 10 m and in the X direction is 20 m. The soil samples were prepared for analysis, and some physical and chemical analyses were made in the samples by routine methods. The data concerning with soil P levels were analyzed comparatively according to GIS based interpolation methods of Ordinary Kriging (OK), Simple Kriging (SK) and Universal Kriging (UK). The interpolation methods were also tested with varied semivariogram models of Spherical, Exponential and Gaussian. As a result of cross validations, the best optimal method was found to be interpolation method of UK (Universal RMSE, +/- 0.472) with semivariogram model of guassian for topsoil, whereas it was found to be interpolation method of SK (Simple RMSE, +/-0.323) with semivariogram model of exponential for topsoil. Predicted P values were significantly (p< 0.01) correlated with measured values for both topsoil (r = 0.993) and subsoil (r = 0.980), respectively. Soil P distribution maps were adequately performed by using selected kriging interpolation methods and suitable sernivariogram models. The results indicated that monitoring of site specific P variability on the apple orchard using these GIS based interpolation methods will help to create the effective schemes for agricultural chemical managements such as P fertilization resulting in optimal yield and quality with reduced environmental pollution.