The Spiritwood Aquifer system in North Dakota is conceptualized as a hierarchical arrangement of lithofacies. The indicator geostatistics formalism is being used to characterize and model proportion, geometry, and juxtaposition of sand and gravel (s) and mud and diamicton (m) lithofacies within assemblages that have high proportion of s (aquifers) and low proportion of s (background). At the larger scale (regional complex of facies assemblages), the proportion (P), geometry, and juxtaposition of aquifer assemblages within the background material are also statistically characterized and modeled in the same way. In results thus far, the mean P_m in aquifer assemblages is 0.15 (volume fraction) and the standard deviation is 0.18 as it varies vertically. The thickness of facies m (l_m) in aquifer assemblages has a mean of 6.63 m and standard deviation of 5.12 m and the l_s has a mean of 9.93 m and a standard deviation of 6.3 m. Vertical auto-transition probability models are spherical as developed with these statistics. Lateral correlation indicates a mean lateral extent of less than 300 m for facies m and 1200 m for facies s, with significant variation. The background assemblage has a mean P_m of 0.90 and a standard deviation of 0.119. The mean l_m is 16.93 m with a standard deviation of 20.16 m, and the mean l_s is 3.81 m with standard deviation of 4.86 m. The corresponding auto-transition probability models are exponential. At the larger scale the mean proportion of aquifer assemblages is 0.83 with standard deviation of 0.375. The mean thickness of aquifer assemblages is 14.72 m with a standard deviation of 10.70 m and the mean thickness of the background assemblages is 29.78 m with standard deviation of 28.47 m. Transition probability models for the assemblages at the scale of the regional complex are currently being developed.