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?
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.