We revise a method by Kalnins, Kress and Miller (2010) for constructing a canonical form for symmetry operators of arbitrary order for the Schrodinger eigenvalue equation H psi equivalent to (Delta(2) + V)psi = E psi on any 2D Riemannian manifold, real or complex, that admits a separation of variables in some orthogonal coordinate system. We apply the method, as an example, to revisit the Tremblay and Winternitz (2010) derivation of the Painleve VI potential for a 3rd order superintegrable flat space system that separates in polar coordinates and, as new results, we give a listing of the possible potentials on the two-sphere that separate in spherical coordinates and all two-hyperbolic (two-sheet) potentials separating in horocyclic coordinates. In particular, we show that the Painleve VI potential also appears for a 3rd order superintegrable system on the two-sphere that separates in spherical coordinates, as well as a 3rd order superintegrable system on the two-hyperboloid that separates in spherical coordinates and one that separates in horocyclic coordinates. Our aim is to develop tools for analysis and classification of higher order superintegrable systems on any 2D Riemannian space, not just Euclidean space.