During the last decade one has witnessed an increasing interest in assessing health risks caused by exposure to contaminants present in the soil, air, and water. A key component of any exposure study is a reliable model for the space-time distribution of pollutants. This paper compares the performances of multi-Gaussian and indicator kriging for modeling probabilistically the spatial distribution of arsenic concentrations in groundwater of southeast Michigan, accounting for arsenic data collected at private residential wells and the hydrogeochemistry of the area. The arsenic data set, which was provided by the Michigan Department of Environmental Quality (MDEQ), includes measurements collected between 1993 and 2002 at 8212 different wells. Factorial kriging was used to filter the short-range spatial variability in arsenic concentration, leading to a significant increase (17--65%) in the proportion of variance explained by secondary information, such as type of unconsolidated deposits and proximity to Marshall Sandstone subcrop. Cross validation of well data shows that accounting for this regional background does not improve the local prediction of arsenic, which reveals the presence of unexplained sources of variability and the importance of modeling the uncertainty attached to these predictions. Slightly more precise models of uncertainty were obtained using indicator kriging. Well data collected in 2004 were compared to the prediction model and best results were found for soft indicator kriging which has a mean absolute error of 5.6 ¿g/L. Although this error is large with respect to the USEPA standard of 10 ¿g/L, it is smaller than the average difference (12.53 ¿g/L) between data collected at the same well and day, as reported in the MDEQ data set. Thus the uncertainty attached to the sampled values themselves, which arises from laboratory errors and lack of information regarding the sample origin, contributes to the poor accuracy of the geostatistical predictions in southeast Michigan.