The complex distribution characteristics of soil textures at a large or regional scale are difficult to understand with the current state of knowledge and limited soil profile data. In this study, an indicator variogram was used to describe the spatial structural characteristics of soil textures of 139 soil profiles. The profiles were 2 m deep with sampling intervals of 0.05 m, from an area of 15 km² in the North China Plain. The ratios of nugget-to-sill values (SH) of experimental variograms of the soil profiles in the vertical direction were equal to 0, showing strong spatial auto-correlation. In contrast, SH ratios of 0.48–0.81 in the horizontal direction, with sampling distances of ~300 m, showed weaker spatial auto-correlation. Sequential indicator simulation (SIS) and indicator kriging (IK) methods were then used to simulate and estimate the 3D spatial distribution of soil textures. The outcomes of the 2 methods were evaluated by the reproduction of the histogram and variogram, and by mean absolute error of predictions. Simulated results conducted on dense and sparse datasets showed that when denser sample data are used, complex patterns of soil textures can be captured and simulated realisations can reproduce variograms with reasonable fluctuations. When data are sparse, a general pattern of major soil textures still can be captured, with minor textures being poorly simulated or estimated. The results also showed that when data are sufficient, the reproduction of the histogram and variogram by SIS was significantly better than by the IK method for the predominant texture (clay). However, when data are sparse, there is little difference between the 2 methods.