Summary: We study numerically the performance properties of a class of approximation schemes for systems of convection-diffusion-reaction models with small diffusion. A coupling of the equations by first order and zero order terms is admitted. Higher order conforming finite element methods are applied to minimize the effects of numerical diffusion and artificial mixing of species. To reduce spurious oscillations close to sharp layers or interfaces, streamline upwind Petrov-Galerkin stabilization and shock capturing as an additional stabilization in the crosswind direction are used. In applications of practical interest the reliability and accuracy of the approach is demonstrated.