Tracer methods are widespread in computational geodynamics for modeling the advection of chemical data. However, they present certain numerical challenges, especially when used over long periods of simulation time. We address two of these in this work: the necessity for mass conservation of chemical composition fields and the need for the velocity field to be pointwise divergence free to avoid gaps in tracer coverage. We do this by implementing the hybrid discontinuous Galerkin (HDG) finite element method combined with a mass-conserving constrained projection of the tracer data. To demonstrate the efficacy of this system we compare it to other common finite element formulations of the Stokes system and projections of the chemical composition. We provide a reference of the numerical properties and error convergence rates which should be observed by using these various discretization schemes. This serves as a tool for verification of existing or new implementations. We summarize these data in a reproduction of a published Rayleigh–Taylor instability benchmark, demonstrating the importance of careful choices of appropriate and compatible discretization methods for all aspects of geodynamics simulations.