Abstract: A robust and efficient simulation technique is developed based on the extension of the mimetic finite difference method to multiscale hierarchical hexahedral (corner-point) grids via use of the multiscale mixed finite element method. The implementation of the mimetic subgrid discretization method is compact and generic for a large class of grids, and thereby, suitable for discretizations of reservoir models with complex geologic architecture. Flow equations are solved on a coarse grid where basis functions with subgrid resolution accurately account for subscale variations from an underlying fine-scale geomodel. The method relies on the construction of approximate velocity spaces that are adaptive to the local properties of the differential operator. A variant of the method for computing velocity basis functions is developed that utilizes an adaptive local-global algorithm to compute multiscale velocity basis functions by capturing the principal characteristics of global flow. Both local and local-global methods generate subgrid-scale velocity fields that reproduce the impact of fine-scale stratigraphic architecture. By using multiscale basis functions to discretize the flow equations on a coarse grid, one can retain the efficiency of an upscaling method, while at the same time produce detailed and conservative velocity fields on the underlying fine grid. A variant of the multiscale method is also implemented for modeling subsurface reactive flow. The accuracy and efficacy of the multiscale method is compared to that of fine-scale models and of coarse-scale models with no subgrid treatment for several subsurface fluid-flow scenarios. Numerical experiments involving two-phase flow and transport phenomena are carried out on high-resolution corner-point grids that explicitly represent example stratigraphic architectures found in real-life shallow marine and turbidite reservoirs. The multiscale method is several times faster than the direct solution of the fine-scale problem and yields more accurate solutions than coarse-scale modeling techniques that resort to explicit effective properties. The accuracy of the multiscale simulation method with adaptive local-global velocity basis functions are compared to that of the local velocity basis functions. The multiscale simulation results are consistently more accurate when the local-global method is employed for computing the velocity basis functions.