If $B(p, r)$ is a small geodesic ball, then $\mathrm{Vol}(\overline{B}(p, r)) = \int_0^r \mathrm{Area}(\partial B(p, \varepsilon)) d\varepsilon$ 
Let $(M, g)$ be an $n$-dimensional Riemannian manifold, and let $\delta > 0$ be such that $B(p, \delta)$ is a geodesic ball around $p$, i.e. it is equal to $\mathrm{exp}_p(B(0, \delta))$. On this ball, we can define the Riemannian volume form $\mathrm{vol}_g$, and compute volumes of smaller closed balls.


Prove that, for $r \in (0, \delta)$, we have $$\mathrm{Vol}(\overline{B}(p, r)) = \int_0^r \mathrm{Area}(\partial B(p, \varepsilon)) d\varepsilon, $$ where the left hand side is the volume of the ball, and $\mathrm{Area}(\partial B(p, \varepsilon))$ is the $(n-1)$ dimensional volume of this geodesic sphere, computed using the induced Riemannian volume form (i.e. the volume form $i_N \mathrm{vol}_g$, where $N$ is the ourward pointing unit normal vector field to this submanifold and $i_N$ is interior multiplication).

By diffeomorphism invariance, we have that $$\mathrm{Vol}(\overline{B}(p, r)) = \int_{\overline{B}(p, r)} \mathrm{vol}_g = \int_{\overline{B}(0, r)} \left( \mathrm{exp}_p \right)^*(\mathrm{vol}_g), $$ where the last integral is inside $T_pM$.
By the Gauss lemma, we also have that the outward pointing unit normal vector field to small geodesic spheres is $$\partial_r = \nabla r = \frac{x^i}{r} \partial_i, $$ where $(x^1, \cdots, x^n)$ are the Riemannian normal coordinates in this geodesic ball, and $$r(x^1, \cdots, x^n) = \sqrt{(x^1)^2 + \cdots + (x^n)^2} $$ is the radial distance function (which is equal to the Riemannian distance from $p$). Hence, $$ \mathrm{Area}(\partial B(p, \varepsilon)) = \int_{\partial B(p, \varepsilon)} i_{\partial_r} \mathrm{vol}_g = \int_{\partial B(0, \varepsilon)} \left( \mathrm{exp}_p \right)^*(i_{\partial_r} \mathrm{vol}_g), $$ where, again, the last integral is inside $T_pM$.
From this question (Interior Product and Pullback Properties), we know that $$\left( \mathrm{exp}_p \right)^*(i_{\partial_r} \mathrm{vol}_g) = i_{\left( \mathrm{exp}_p \right)^* \partial_r} (\left( \mathrm{exp}_p \right)^* \mathrm{vol}_g). $$
As such, we need to show that $$\int_{\overline{B}(0, r)} \left( \mathrm{exp}_p \right)^*(\mathrm{vol}_g) = \int_0^r \int_{\partial B(0, \varepsilon)} i_{\left( \mathrm{exp}_p \right)^* \partial_r} (\left( \mathrm{exp}_p \right)^* \mathrm{vol}_g) d\varepsilon, $$ but I am not sure how to prove this. Should I try to expand all of these volume forms from above into their coordinates representations?
 A: You can further "decompose" your integrals using polar coordinates (these are known are Riemannian polar coordinates). More precisely, let  $$\Theta : [0, \pi]^{n-2} \times [0, 2\pi] \times [0, r) \to B(0, r) $$ be the polar coordinate chart of $\mathbb{R}^n$. Explicitly, this is written as $$\Theta (\theta_1, \cdots, \theta_{n-2}, \theta_{n-1}, s) =
\begin{pmatrix}
s\cos(\theta_1) \\ s\sin(\theta_1) \cos(\theta_2) \\ \vdots \\ s\sin(\theta_1)\cdots \sin(\theta_{n-2})\cos(\theta_{n-1}) \\ 
s\sin(\theta_1)\cdots \sin(\theta_{n-2})\sin(\theta_{n-1})
\end{pmatrix}. $$
One can then show (by just applying the change of variables formula and Fubini's theorem), that $$\int_{B(0, r)} f(x)dx = \int_0^r d\varepsilon \int_{\partial B(0, \varepsilon)} f(s) d\mathrm{vol}_{\partial B(0, \varepsilon)}(s), $$ where $d\mathrm{vol}_{\partial B(0, \varepsilon)}$ is the induced Riemannian volume form on the sphere, or, written explicitly, $$d\mathrm{vol}_{\partial B(0, \varepsilon)} = i_{\partial_r} (dx^1 \wedge \cdots \wedge dx^n). $$
Let us now use the above fact to answer your question.
In Riemannian normal coordinates, we have $$\left(\mathrm{exp}_p \right)^* \mathrm{vol}_g = \sqrt{\det(g_{ij})} dx^1 \wedge \cdots \wedge dx^{n}, $$ where $g = g_{ij} dx^i dx^j$ is the metric in these coordinates.
Hence, applying the above identity, we have $$\mathrm{Vol}(B(0, r)) = \int_0^r d\varepsilon \int_{\partial B(0, \varepsilon)} \sqrt{\det(g_{ij})} i_{\partial_r} (dx^1 \wedge \cdots \wedge dx^n). $$
This takes care of the volume of the ball. Now, let us proceed similarly for the area of the sphere. As you have shown,
$$ \mathrm{Area}(\partial B(p, \varepsilon)) = \int_{\partial B(0, \varepsilon)} \left( \mathrm{exp}_p \right)^*(i_{\partial_r} \mathrm{vol}_g) = \int_{\partial B(0, \varepsilon)} i_{\left( \mathrm{exp}_p \right)^* \partial_r} (\left( \mathrm{exp}_p \right)^* \mathrm{vol}_g). $$
Since interior multiplication is $C^\infty(M)$ linear, we have $$\int_{\partial B(0, \varepsilon)} i_{\left( \mathrm{exp}_p \right)^* \partial_r} (\left( \mathrm{exp}_p \right)^* \mathrm{vol}_g) = \int_{\partial B(0, \varepsilon)} \sqrt{\det(g_{ij})} i_{\left( \mathrm{exp}_p \right)^* \partial_r} (dx^1 \wedge \cdots \wedge dx^n). $$
Lastly, we have that $$ \left(\mathrm{exp}_p \right)^* \partial_r = \partial_r, $$ since, in Riemannian normal coordinates, the exponential map reduces to the identity map (or, more precisely, writing $\partial_r$ is coordinates is the same as computing $ \left( \mathrm{exp}_p^{-1} \right)^* (\partial_r)$ with this second $\partial_r$ being the usual radial derivative in $\mathbb{R}^n$). This concludes the proof.
