The development of numerical algorithms for modeling flow processes in large-scale highly heterogeneous formations is very challenging because the properties of natural geologic porous formations (e.g., permeability) display high variability levels and complex spatial correlation structures, which span a rich hierarchy of length scales. Thus, it is usually necessary to resolve a wide range of length and time scales in order to obtain accurate predictions of the flow, mechanical deformation, and transport processes under investigation. In practice, however, some type of coarsening (or upscaling) of the detailed model is usually performed before the model can be used to simulate complex displacement processes. Many approaches have been developed and applied successfully when a scale separation adequately describes the spatial variability of the subsurface properties (e.g., permeability) that have bounded variations. The quality of these approaches deteriorates for complex heterogeneities, especially when the contrast in the media properties is large, e.g., in the case of fractured porous media. In this talk, I will describe multiscale model reduction techniques that can be used in upscaling flow equations as well as in domain decomposition methods when media properties have high contrast and are spatially heterogeneous. Numerical results will be presented that show that one can improve the accuracy of multiscale methods and obtain contrast-independent preconditioners.