Introduction MotivationBuoyancy is the primary driving force behind convection in the Earth's mantle. Contributing factors to buoyancy in the mantle include lateral contrasts in temperature and composition. In the case of thermochemical flows, mantle buoyancy depends upon both of these factors. Heat sources and sinks affecting the thermal state of the Earth's mantle include radiogenic heating, heating from the outer core, and cooling at the Earth's surface (Turcotte & Schubert, 2002). Thermally driven buoyancy instabilities can arise when the rate of thermal transport via advection exceeds that of diffusion. Such situations include the rise of hot upwellings from the core and the descent of cold downwellings from the surface (e.g., a subducting slab). Buoyancy instabilities can also be caused by lateral variations in thermal boundary conditions. Lateral contrasts in mantle composition can occur for several reasons, including transitions between oceanic and continental lithosphere, rapid subduction of oceanic lithosphere, and deep dense compositional piles. The Earth's Large Low Shear wave Velocity Provinces may also be influenced by thermal and compositional gradients (Davies et al., 2015;McNamara, 2019).Geophysical flows involving sharp compositional contrasts are notoriously difficult to model numerically. Challenges include both spurious oscillations and extraneous diffusion (Lenardic & Kaula, 1993). Numerical methods employed to minimize these errors include the use of particles (Tackley & King, 2003), level sets (Hillebrand et al., 2014), and hybrid methods (Samuel & Evonuk, 2010).Another even more fundamental challenge is to ensure the software that implements the numerical solution has been coded correctly. By definition, model developers produce code to solve problems for which solutions are unknown. The verification and validation process in modeling and simulation often involves qualitative

In this study, we manufacture an exact solution for a set of 2D thermochemical mantle convection problems. The derivation begins with the specification of a stream function corresponding to a non-stationary velocity field. The method of characteristics is then applied to determine an expression for composition consistent with the velocity field. The stream function formulation of the Navier-Stokes equation is then applied to solve for temperature. The derivation concludes with the application of the advection-diffusion equation for temperature to solve for the internal heating rate consistent with the velocity, composition, and temperature solutions. Due to the large number of terms, the internal heating rate is computed using Maple which is also made available in Fortran. Using the method of characteristics allows the compositional transport equation to be solved without the addition of diffusion or source terms. As a result, compositional interfaces remain sharp throughout time and space in the exact solution. The exact solution presented allows for precision testing of thermochemical convection codes for correctness and accuracy.