The characterization of the soil habitat is of fundamental importance to an understanding of processes associated with sustainable management such as environmental flows, bioavailability, and soil ecology. We describe a method for quantifying and explicitly modeling the heterogeneity of soil using a stochastic approach. The overall aim is to develop a model capable of simultaneously reproducing the spatial statistical properties of both the physical and biological components of soil architecture. A Markov chain Monte Carlo (MCMC) methodology is developed that uses a novel neighborhood and scanning scheme to model the two-dimensional spatial structure of soil, based on direct measurements made from soil thin sections. The model is considerably more efficient and faster to implement than previous approaches, and allows accurate modeling of larger structures than has previously been possible. This increased efficiency also makes it feasible to extend the approach to three dimensions and to simultaneously study the spatial distribution of a greater number of soil components. Examples of two-dimensional structures created by the models are presented and their statistical properties are shown not to differ significantly from those of the original visualizations.