We describe a multilevel multiscale mimetic (M3) method for solving two-phase flow (water and oil) in a heterogeneous reservoir. The governing equations are the elliptic equation for the reservoir pressure and the hyperbolic equation for the water saturation. On each time step, we first solve the pressure equation and then use the computed flux in an explicit upwind finite volume method to update the saturation. To reduce the computational cost, the pressure equation is solved on a much coarser grid than the saturation equation. The coarse-grid pressure discretization captures the influence of multiple scales via the subgrid modeling technique for single-phase flow recently proposed in [?]. We extend significantly the applicability of this technique by developing a new robust and efficient method for estimating the flux coarsening parameters. Specifically, with this advance the M3 method can handle full permeability tensors and general coarsening strategies, which may generate polygonal meshes on the coarse grid. These problem dependent coarsening parameters also play a critical role in the interpolation of the flux, and hence, in the advection of saturation for two-phase flow. Numerical experiments for two-phase flow in highly heterogeneous permeability fields, including layer 68 of the SPE Tenth Comparative Solution Project, demonstrate that the M3 method retains good accuracy for high coarsening factors in both directions, up to 64 for the considered models. Moreover, we demonstrate that with a simple and efficient temporal updating strategy for the coarsening parameters, we achieve accuracy comparable to the fine-scale solution, but at a fraction of the cost.