Since the 1950s, numerical models are widely used in reservoir simulators to predict and optimize oil recovery from petroleum reservoirs. Commercial simulators typically combine multiphase porous media flow models with a separate module for the transport of tracers. By decoupling flow and transport, the transport equation can be solved using a more accurate and efficient numerical scheme than the fully implicit first order scheme that is commonly used to solve the flow equations. The accuracy of the numerical transport scheme is of high importance in Enhanced Oil Recovery modelling in order to accurately predict the influence of EOR techniques on the oil production. Of all numerical methods considered in this research, the explicit high-resolution flux-limiter method with the van Leer limiter performed best in terms of accuracy and efficiency. The accuracy and monotonicity of the numerical transport scheme strongly depend on the underlying flowsolution. To ensure monotonicity of the scheme in all model situations considered, a partially implicit method is introduced that switches to the monotone first order implicit scheme where necessary. In addition, care must be taken when modelling influences of the polymer and surfactant concentration on the flow. Bad choices can lead to instabilities or cause deviations in the model outcome. The conventional approach for modelling the hydrodynamic acceleration of polymer via a constant velocity enhancement factor results in an ill-posed system. To obtain a well-posed system, some form of a saturation dependent velocity enhancement factor could be used instead. However, which physical and numerical model is best suited for modelling the hydrodynamic acceleration remains an open question. When modelling surfactant flooding, the implementation of a discontinuous transition in relative permeabilities around the surfactant front results in oscillations in the numerical solution over time. These oscillations can be resolved or at least diminished by applying some form of interpolation between the two sets of relative permeabilities. However, there is no guarantee that the oscillations disappear completely and the interpolation function has to be carefully selected as it can have a large impact on the solution profile.