We present a global/local model reduction for fast multiscale reservoir simulations in highly heterogeneous porous media. Our approach identifies a low-dimensional structure in the solution space. We introduce an auxiliary variable (the velocity field) in our model reduction that achieves a high compression of the model. This compression is achieved because the velocity field is conservative for any low-order reduced model in our framework, whereas a typical global model reduction that is based on proper-orthogonal-decomposition (POD) Galerkin projection cannot guarantee local mass conservation. The lack of mass conservation can be observed in numerical simulations that use finite-volume-based approaches. The discrete empirical interpolation method (DEIM) approximates fine-grid nonlinear functions in Newton iterations. This approach delivers an online computational cost that is independent of the fine-grid dimension. POD snapshots are inexpensively computed with local model-reduction techniques that are based on the generalized multiscale finite-element method (GMsFEM) that provides (1) a hierarchical approximation of the snapshot vectors, (2) adaptive computations with coarse grids, and (3) inexpensive global POD operations in small dimensional spaces on a coarse grid. By balancing the errors of the global and local reduced-order models, our new methodology provides an error bound in simulations. Our numerical results, by use of a two-phase immiscible flow, show a substantial speedup, and we compare our results with the standard POD-DEIM in a finite-volume setup.