P>Accurate predictions of stress distribution in the lithosphere are of
major importance for approaching more realistic numerical models of the
mantle-lithosphere system. Since stress fields in the lithosphere
computed in convection models differ substantially between viscous and
viscoelastic rheologies, it is essential to employ a viscoelastic
rheology when accurate stresses are to be predicted in mantle convection
models involving the lithosphere. This difference in stress distribution
and magnitude has important implications for accurate modelling of
stress-dependent processes like power-law creep, shear heating and
plasticity. A further requirement for computation of accurate stress
fields is to ensure numerically divergence-free solutions in the
Boussinesq approximation. We present the technical background required
for implementation of numerically incompressible solutions and for
implementation of a Maxwell viscoelastic rheology in the frame of the
finite element method (FEM). We employ the Jaumann invariant stress
derivative in our implementation and demonstrate that the choice of the
invariant stress derivative is irrelevant for geodynamic simulations. We
discuss potential numerical advantages of a viscoelastic rheology when
large viscosity variations are applied in thermal convection models. Due
to the physical transition from effectively viscous to elastic behaviour
in a viscoelastic model, the introduction of viscosity cut-offs
generally applied in viscous models can be avoided.