  Motivated by challenges in Earth mantle convection, we present a massively
parallel implementation of an Eulerian-Lagrangian method for the
advection-diffusion equation in the advection-dominated regime. The advection
term is treated by a particle-based, characteristics method coupled to a
block-structured finite-element framework. Its numerical and computational
performance is evaluated in multiple, two- and three-dimensional benchmarks,
including curved geometries, discontinuous solutions, pure advection, and it is
applied to a coupled non-linear system modeling buoyancy-driven convection in
Stokes flow. We demonstrate the parallel performance in a strong and weak
scaling experiment, with scalability to up to $147,456$ parallel processes,
solving for more than $5.2 \times 10^{10}$ (52 billion) degrees of freedom per
time-step.
