Hi Till,
it’s very unlikely that there is a global orthogonal parametrization. Though, it’s straightforward to fix a parameter and to determine two parameters such that they are orthogonal to the fixed parameter but not to each other (Cox & Reid, 1987). One such solution: geometric mean is orthogonal to scale and index parameter. Taking scale as fixed has no analytic solution for replacing the index parameter. I’m currently writing a manuscript with these results.
As you already found out in your simulations the mean is orthogonal to the shape parameter (of the gamma distribution). So you may try to look for the second orthogonal parameter to replace the index parameter.
If you wish to try a parametrization with a geometric mean you will find my implementation of inverse digamma function helpful.
Cox, D. R., & Reid, N. (1987). Parameter orthogonality and approximate conditional inference. Journal of the Royal Statistical Society: Series B (Methodological) , 49 (1), 1-18.