Fatigue lives of rolling element bearings exhibit a wide scatter due to the statistical nature of the rolling contact fatigue failure process. Empirical life models that account for this dispersion do not provide insights into the physical mechanisms that lead to this scatter. One of the primary reasons for dispersion in lives is the stochastic nature of the bearing material. Here, a damage mechanics based fatigue model is introduced in conjunction with the idea of discrete material representation that takes the effect of material microstructure explicitly into account. Two sources of material randomness are considered: (1) the topological randomness due to geometric variability in the material microstructure and (2) the material property randomness due to nonuniform distribution of properties throughout the material. The effect of these variations on the subsurface stress fields in rolling element line contacts is studied. The damage model, which incorporates cyclic damage accumulation and progressive degradation of material properties with rolling contact cycling, is used to study the mechanisms of subsurface initiated spalling in bearing contacts. Crack initiation as well as propagation stages are modeled using damaged material zones in a unified framework. The spalling phenomenon is found to occur through microcrack initiation below the surface where multiple microcracks coalesce and subsequent cracks propagate to the surface. The computed crack trajectories and spall profiles are found to be consistent with experimental observations. The microcrack initiation phase is found to be only a small fraction of the total spalling life and the scatter in total life is primarily governed by the scatter in the propagation phase of the cracks through the microstructure. Spalling lives are found to follow a three-parameter Weibull distribution more closely compared to the conventionally used two-parameter Weibull distribution. The Weibull slopes obtained are within experimentally observed values for bearing steels. Spalling lives are found to follow an inverse power law relationship with respect to the contact pressure with a stress-life exponent of 9.35.