Abstract: Soil pH value is the major property influencing the availability of soil nutrients and behaviors of agricultural chemicals. Understanding the spatial distribution of soil pH values has vital significance for nutrient management and soil pollution prediction. Geostatistical method was used to study the distribution of soil pH values in a watershed in the gully region of the Loess Plateau where landforms, landuse patterns and soil types vary greatly. The results showed that soil pH values in the watershed presented the Spherical-Exponential combined variogram models. The spatial heterogeneity was mainly derived from the random factors such as landuse patterns and soil erosion in the watershed. Crossed Kriging on soil pH value, using organic matter as auxiliary variables, was consistently more effective than Kriging on soil pH value alone and could be used to estimate soil pH value distribution in watershed. The error of crossed Kriging was existed in the complicated environmental factors. Simulated values by sequential Gaussian conditional simulation had the same range but lower average soil pH values compared with measured data. The error was derived from the special Kriging method and the Gaussian treat during the simulation.