Summary: This paper presents higher order methods for the numerical modeling of two-phase flow with simultaneous transport and adsorption of viscosifying species within the individual phases in permeable porous media. The numerical scheme presented addresses the three major challenges in simulating this process. Firstly, the component transport is strongly coupled with the viscous and capillary forces that act on the movement of the carrier phase. The discretization of the capillary parts is especially difficult since its effect on flow yields non-linear parabolic conservation equations. These are amenable to non-linear finite elements (FEs), while the capillary contribution on the component transport is first-order hyperbolic, where classical FEs are unsuitable. We solve this efficiently by a Strang splitting that uses finite volumes (FVs) with explicit time-stepping for the viscous parts and a combined finite element-finite volume (FEFV) scheme with implicit time-stepping for the capillary parts. Secondly, the components undergo hydrodynamic dispersion and discerning between numerical and physical dispersion is essential. We develop higher-order formulations for the phase and component fluxes that keep numerical dispersion low and combine them with implicit FEs such that the non-linearities of the dispersion tensor are fully incorporated. Thirdly, subsurface permeable media show strong spatial heterogeneity, with coefficients varying over many orders of magnitude and geometric complexity that make the use of unstructured grids essential. In this work, we employ node-centered FVs that combine their ability to resolve flow with the flexibility of FEs. Numerical examples of increasing complexity are presented that demonstrate the convergence and robustness of our approach and prove its versatility for highly heterogeneous, and geometrically complex fractured porous media.