Fubini's formula on fibered Riemannian manifolds. Let $\pi:(M,g)\to(N,h)$ be a Riemannian submersion (that is, $\pi$ is a submersion satisfying $g|_{\ker(d\pi)^\perp}=\pi^*h|_{\ker(d\pi)^\perp}$). Each fiber of $\pi$ is a Riemannian submanifold of $M$, given by $(F_y,g_y):=\left(\pi^{-1}(y),g|_{\pi^{-1}(y)}\right)$. In the case that $\pi$ is a canonical projection on a Riemannian product manifold, and $f\in C^\infty_cM$ is a compactly supported smooth function, we can write Fubini's formula suggestively as
$$
\int_Mf(x)d\mu_g(x)=\int_N\left(\int_{F_y}f(x)d\mu_{g_y}(x)\right)d\mu_h(y)
$$
Where $\mu_g$ is the measure on $M$ induced by $g$, etc. Intuitively, it seems that this formula also should hold for a general Riemannian submersion $\pi$, but it is not so clear how to show this. $M$ need not look like a Riemannian product space, even locally, since $\ker(d\pi)^\perp$ need not be integrable. As a result, $g$ does not decompose nicely in adapted coordinates, making a direct computation difficult. I suspect there's a more straightforward geometrical argument, or a simple counterexample that I'm missing.
 A: It seems that the missing piece was was the use of an adapted orthonormal coframe (thanks to Ted Shifrin's suggestion) and a version of Fubini's theorem for differential forms, described in this answer. Let $\nu$ be an $n$-form on $N$, and $\mu$ be an $(m-n)$-form on $M$ (with both smooth and $\mu$ compactly supported). With suitable choices of orientation, we have
$$
\int_M\pi^*\nu\wedge\mu=\int_N\left(\int_{\pi^{-1}(y)}\mu\right)\nu(y)
$$
This can be shown straightforwardly by working on adapted local coordinates (and the hypotheses can be weakened significantly, but that isn't needed here).
We can w.l.o.g. consider the case where $\operatorname{supp}(f)$ is contained in the domain $U\subseteq M$ of a local adapted orthonormal frame $e^1,\cdots,e^m$, where $e^i=\pi^*(\widetilde{e}^i)$, and $\widetilde{e}^1,\cdots,\widetilde{e}^n$ are a local orthonormal frame on $\pi(U)\subseteq N$. Using the shorthand $e^{k\cdots l}:=e^k\wedge\cdots\wedge e^l$, and choosing appropriate orientations, we can write the integral in terms of these frames.
\begin{align*}
\int_Uf(x)d\mu_g(x)=&\int_Uf(x)e^{1\cdots m}(x) \\
=&\int_U\pi^*\left(\widetilde{e}^{1\cdots n}(x)\right)\wedge\left(f(x)e^{(n+1)\cdots m}(x)\right)
\end{align*}
This expression can be Fubinated, yielding the desired result.
\begin{align*}
=&\int_{\pi(U)}\left(\int_{U\cap F_y}f(x)e^{(n+1)\cdots n}(x)\right)\widetilde{e}^{1\cdots n}(y) \\
=&\int_{N}\left(\int_{F_y}f(x)d\mu_{g_y}(x)\right)d\mu_h(y) \\
\end{align*}
