A novel method for simulating multi-phase flow in porous media is presented. The approach is based on a
control volume finite element mixed formulation and new force-balanced finite element pairs. The novelty of
the method lies in: (a) permitting both continuous and discontinuous description of pressure and saturation
between elements; (b) the use of arbitrarily high-order polynomial representation for pressure and velocity
and (c) the use of high-order flux-limited methods in space and to time avoid introducing non-physical
oscillations while achieving high-order accuracy where and when possible. The model is initially validated
for two-phase flow. Results are in good agreement with analytically obtained solutions and experimental
results. The potential of this method is demonstrated by simulating flow in a realistic geometry composed of
highly permeable meandering channels.