Even on modern supercomputer architectures, Earth mantle simulations are so compute intensive that they are considered grand challenge applications. The dominating roadblocks in this branch of Geophysics are model complexity and uncertainty in parameters and data, e.g., rheology and seismically imaged mantle heterogeneity, as well as the enormous space and time scales that must be resolved in the computational models. This article reports on a massively parallel all-at-once multigrid solver for the Stokes system as it arises in mantle convection models. The solver employs the hierarchical hybrid grids framework and demonstrates that a system with coupled velocity components and with more than a trillion ($1.7 \cdot 10^{12}$) degrees of freedom can be solved in about 1,000 s using 40,960 compute cores of JUQUEEN. The simulation framework is used to investigate the influence of asthenosphere thickness and viscosity on upper mantle velocities in a static scenario. Additionally, results for a time-dependent simulation with a time-variable temperature-dependent viscosity model are presented.