Reservoirs are inherently heterogeneous. The depositional environment, compaction, deformation and cementation can lead to spatial variability in petrophysical properties, over many orders of magnitude and length-scales, with a profound effect on the flow behaviour of the fluids flowing within it. Conventional reservoir modelling workflows have adopted k-orthogonal meshes that are fixed in time and space. Due to the complex geometries found in reservoirs and the multi-scale structure of the governing equations, a more flexible approach is needed. This thesis presents an alternative to conventional grid-based modelling and simulation workflows. SBM is used to capture geological heterogeneity at all scales that impact flow, without reference to an underlying grid. Petrophysical properties within these domains are modelled as internally homogeneous. We show that cell-to-cell variability in reservoir models is not needed to capture multi-phase flow. Instead one can use a small number of larger, internally homogeneous but geometrically complex geological domains. For flow simulations, geomodels are discretised using unstructured tetrahedral meshes and the governing equations are solved using a CVFEM. The proposed simulator incorporates DAMO, designed to tackle the unique challenges of spatially varying petrophysical properties in porous media. Multi-scale features are better captured at lower computational cost by refining the mesh resolution where/when necessary and coarsening elsewhere. Material discontinuities in porous media strongly influence multi-phase flow. Traditional CVFEM smear out jump discontinuities of saturation or concentration. The high aspect ratio elements typically required to discretise the heterogeneous domains in subsurface reservoirs are also very challenging on standard CVFEM. To overcome these issues we develop a new discretisation based on the bubble enriched FEM. This new method is able to accurately represent saturation discontinuities on geologically complex models.