The spatial distribution information of clay layer occurrence depth (CLOD), particularly the spatial distribution maps of occurrence of clay layers at depths less than a certain threshold, in alluvial soils is crucial to designing appropriate plans and measures for precision agriculture and environmental management in alluvial plains. Markov chain geostatistics (MCG), which was proposed recently for simulating categorical spatial variables, can objectively decrease spatial uncertainty and consequently increase prediction accuracy in simulated results by using nonlinear estimators and incorporating various interclass relationships. In this paper, a MCG method was suggested to simulate the CLOD in a meso-scale alluvial soil area by encoding the continuous variable with several threshold values into binary variables (for single thresholds) or a multi-class variable (for all thresholds being considered together). Related optimal prediction maps, realization maps, and occurrence probability maps for all of these indicator-coded variables were generated. The simulated results displayed the spatial distribution characteristics of CLOD within different soil depths in the study area, which are not only helpful to understanding the spatial heterogeneity of clay layers in alluvial soils but also providing valuable quantitative information for precision agricultural management and environmental study. The study indicated that MCG could be a powerful method for simulating discretized continuous spatial variables. Copyright © 2009 John Wiley & Sons, Ltd.