Imaging subsurface structures, such as salt domes, magma reservoirs, or subducting plates, is a major challenge in geophysics. Seismic imaging methods are, so far, the most precise methods to open a window into the Earth. However, the methods may not yield the exact depth or size of the imaged feature and may become distorted by phenomena such as seismic anisotropy, fluid flow, or compositional variations. A useful complementary method is therefore to simulate the mechanical behavior of rocks on large timescales, and compare model predictions with observations. Recent studies have used the (nonlinear) Stokes equations and geometries from seismic studies in combination with an adjoint-based approach to invert for rheological parameters that are consistent with surface observations such as GPS velocities. Nevertheless, it would be useful to employ other surface observations, such as principal stress directions, as constraints as well. Here, we derive the adjoint formulation for the case that principal stress directions are used as observables with respect to rheological parameters. Both an algebraic and a discretized derivation of the adjoint equations are described. This thus enables the usage of two data fields - surface velocities and stress directions - as a misfit for the inversion. We test the performance of the inversion for principal stress directions on simplified 3D test cases. Finally, we demonstrate how the adjoint approach can be employed to compute 3D geodynamic sensitivity kernels, which highlight the areas in the model domain that have the largest impact on the misfit value of a particular point. This provides a simple, yet powerful, way to visualize which parts of the model domain are of key importance if changing rheological constants.