The behavior of particulate media, such as sands, is encoded at the granular-scale and hence methods for upscaling such behavior across relevant scales of interest—from granular-scale (~1 mm) to field-scale (>1m)—are needed to attain a more accurate prediction of soil behavior. Multi-scale analysis is especially important under extreme conditions, such as strain localization, penetration, or liquefaction, where the classical constitutive description may no longer apply. In this paper, internally consistent probabilistic models for undrained shear strength and Young's modulus are developed at multiple scales, and incorporated into a simulation framework where refinement of the material description to finer scales is pursued only as necessary. This probabilistic simulation approach is then coupled with the finite element method. Numerical examples are presented to show how the performance of the geosystem is influenced by taking into account multi-scale random fields.