We apply the adjoint method to the post-glacial sea level problem in order to calculate the derivatives of measurements with respect to mantle viscosity and ice sheet history. These derivatives, or kernels, quantify the linearised sensitivity of such measurements to their underlying model parameters. The adjoint method, which has previously been used within a range of other geophysical applications, enables efficient calculation of theoretically exact sensitivity kernels within laterally heterogeneous earth models having a range of linear or non-linear viscoelastic rheologies. The kernels have a number of applications within the inverse problem. Firstly, they can be used within a gradient-based optimisation method to find a model which minimises some data misfit function. The kernels can also be used to quantify the uncertainty in such a model and hence to provide understanding of which parts of the model are well constrained. Finally, they enable construction of measurements which provide sensitivity to a particular part of the model space. We illustrate potential applications of this method for studies of glacial isostatic adjustment through a range of numerical calculations relative to a spherically symmetric background model.