Thermal convection in the Earth's mantle is necessary for explaining a multitude of geophysical and geological processes which we can observe and measure at the surface of the Earth, such as continental drift, perturbations in Earth's gravitational field, and variations in oceanic bathymetry and continental elevation. Over the past decade numerous analyses of these convection-related observables have been carried out in the context of seismic tomography based mantle flow models. One of the major objectives of this modelling has been to infer the rheological structure of the mantle. There exist a multitude of numerical methods for the solutions of planetary scale convection. I will present an overview of the variational formulation for a buoyancy induced mantle flow problem and its solution in spherical geometry using generalized spherical harmonics. This approach yields a very large dense system of equations for desirable spatial resolutions at the surface and thus some numerical and practical considerations will also be discussed.