Two hierarchical levels are used to characterize heterogeneity in the Mahomet aquifer system in East-Central Illinois, towards developing ground water models. The relevant Quaternary stratigraphy includes the Mason and Wedron Groups, the Glasford Formation, and the Banner Formation which includes the Mahomet or Sankoty Sand Member. While these units are dominated by either sand and gravel (s) facies or mud and diamicton (m) facies, none are homogeneous, but are heterogeneous mixtures of these facies. The indicator geostatistics formalism is used to both characterize and model the proportions, geometry, and distribution of units at two scales. At the smaller scale, proportions, geometry, and spatial distribution of lithofacies m and s are characterized within facies assemblages, and if stationary, are modeled with transition probabilities. We have developed such models for facies assemblages representing each of the stratigraphic groups/formations mentioned above. At the larger scale, we will in the same way characterize and model the proportions, geometry, and spatial distribution of the assemblages in the overall complex of sediments. We are currently exploring where the boundaries of facies assemblages lie based upon changes in facies proportions, and comparing these to accepted formation boundaries. In some cases the proportions of facies m and s do not change abruptly at elevations thought to correspond to formation boundaries, leading us to explore the topography of the boundaries of facies assemblages and to question if additional, intermediate facies assemblages need to be defined between the facies assemblages that were expected a-priori.