  The geodesic Markov chain Monte Carlo method and its variants enable
computation of integrals with respect to a posterior supported on a manifold.
However, for regular integrals, the convergence rate of the ergodic average
will be sub-optimal. To fill this gap, this paper extends the efficient
posterior integration method of Oates et al. (2017) to the case of a Riemannian
manifold. In contrast to the original Euclidean case, no non-trivial boundary
conditions are needed for a closed manifold. The method is assessed through
simulation and deployed to compute posterior integrals for an Australian
Mesozoic paleomagnetic pole model, whose parameters are constrained to lie on
the manifold $M = \mathbb{S}^2 \times \mathbb{R}_+$.
