SUMMARYStress concentrations near grain boundaries, precipitates, and similar micro-heterogeneities nucleate instabilities leading to macroscale fracture. As it is not practical to model each flaw explicitly, their ensemble effect is modeled statistically. Accounting for this aleatory uncertainty requires smaller specimens (e.g., small finite elements) to have generally higher and more variable strengths, which is necessary for the initial failure probability of a finite domain to be unaffected by its discretization into elements. Localization itself, which might be attributed to constitutive instability, requires realistic numerical perturbations to predict bifurcations such as radial cracking in axisymmetric problems. These perturbations, stemming from microscale heterogeneity, are incorporated in simulations by imposing statistical spatial variability in the parameters of an otherwise conventional (deterministic and scale-independent) damage model. This approach is attractive for its algorithmic simplicity and straightforward calibration from standard strength tests. In addition, it results in virtually no loss of efficiency or robustness relative to deterministic models and accommodates general three-dimensional loading. Despite these advantages, some significant challenges remain and are discussed. However, it is demonstrated that including aleatory uncertainty with associated scale effects significantly improves predictiveness on large-scale computational domains, where it is impractical to resolve each crack or localization zone.