Numerical modelling of Earth's mantle is a complex, and computationally demanding task due to, amongst others, the broad spectrum of temporal and spatial scales playing a role in mantle flow, large uncertainties in the physical properties of mantle material, with large and localised transitions in viscosity and density. This thesis introduces and analyses a number of numerical techniques that may bring a significant contribution in meeting some of these challenges. Firstly, we introduce a novel time integration scheme for free surface movement in mantle convection models that is more accurate and stable for large time steps. Secondly, we extend the capabilities of anisotropic mesh optimisation, which allows efficient focussing of mesh resolution, to handle cylindrical and spherical shell domains and demonstrate that a significant reduction in the required number of degrees of freedom is possible while maintaing accuracy. Finally, to verify correctness, and evaluate and compare properties of various numerical schemes, we derive an extensive suite of analytical solutions to the Stokes equations governing mantle flow in cylindrical and spherical shell domains, with physically relevant boundary conditions. As a numerical benchmark they also serve to facilitate comparisons of different geodynamical models, and the further development of numerical techniques to improve these.