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?