We develop a local basis model-order reduction technique for approximation of flux/pressure fields based on local proper orthogonal decompositions (PODs) consistently glued together using the multiscale Mixed FEM (MsMFEM) framework on a coarse grid. Based on snapshots from one or more simulation run, we perform SVDs for the flux distribution over coarse grid interfaces and use the singular vectors corresponding the largest singular values as boundary conditions for the multiscale flux basis functions.The span of these basis functions matches (to prescribed accuracy) the span of the snapshots over coarse grid faces. Accordingly, the complementary span (whats left) can be approximated by local PODs on each coarse block giving a second set of local/sparse basis functions. The reduced system unknowns corresponding to the second set of basis functions can be eliminated to keep the system size low. To assess the accuracy, we apply the methodology to a realistic test problem (two-phase compressible flow including gravity) with several wells and compare to results obtained from full order simulations. Both changing well configurations and changing well placements (with local update of bases) are considered. In addition, comparison to standard POD is considered.