We present a new class of numerical methods for solving stochastic differential equations with additive noise on general Riemannian manifolds with high weak order of accuracy. In opposition to the popular approach with projection methods, the proposed methods are intrinsic: they only rely on geometric operations and avoid coordinates and embeddings. We provide a robust and general convergence analysis and an algebraic formalism of exotic planar Butcher series for the computation of order conditions at any high order. To illustrate the methodology, an explicit method of second weak order is introduced, and several numerical experiments confirm the theoretical findings and extend the approach for the sampling of the invariant measure of Riemannian Langevin dynamics.