[proofplan]
We prove the identity first for nonnegative Borel functions inside a single coordinate chart by translating the Riemannian integrals into Euclidean coordinates and applying the [Euclidean Coarea Formula](/page/Coarea%20Formula) for Lipschitz functions. The metric tensor contributes the volume density and the tangential hypersurface Jacobian, and a determinant computation shows that these two factors combine with $|\nabla u|_g$ exactly as required. A countable locally finite [partition of unity](/page/Partition%20of%20Unity) then decomposes the global Borel case into chart-supported pieces. Finally, a Borel-representative argument extends the identity to functions measurable with respect to the completed Riemannian volume measure, and the signed case follows by applying the nonnegative case to the positive and negative parts.
[/proofplan]
[step:Reduce the global Borel identity to chart-supported summands]
We first prove the formula for a nonnegative Borel measurable function $f:M\to[0,\infty]$. Choose a countable smooth atlas $(W_j,\theta_j)_{j\in\mathbb{N}}$ whose coordinate domains cover $M$. Since smooth manifolds are paracompact and second-countable, refine this atlas cover by a countable locally finite open cover $(O_i)_{i\in\mathbb{N}}$ such that each closure $\overline{O_i}$ is contained in some coordinate domain $W_{j(i)}$. The [Smooth Partition of Unity Theorem](/page/Smooth%20Partition%20of%20Unity) gives a countable locally finite smooth partition of unity $(\rho_i)_{i\in\mathbb{N}}$ subordinate to $(O_i)_{i\in\mathbb{N}}$. For each $i$, set $U_i:=W_{j(i)}$ and $\varphi_i:=\theta_{j(i)}|_{U_i}$. Thus $0\leq \rho_i\leq 1$, $\sum_{i=1}^{\infty}\rho_i=1$, and $\operatorname{supp}\rho_i\subset \overline{O_i}\subset U_i$.
For each $i\in\mathbb{N}$, define
\begin{align*}
f_i:M &\to [0,\infty] \\
x &\mapsto \rho_i(x)f(x).
\end{align*}
Since $f_i\geq 0$ and $f=\sum_{i=1}^{\infty}f_i$ pointwise, the [monotone convergence theorem](/theorems/509) gives
\begin{align*}
\int_M f(x)|\nabla u(x)|_g\,d\operatorname{vol}_g(x)
=
\sum_{i=1}^{\infty}
\int_M f_i(x)|\nabla u(x)|_g\,d\operatorname{vol}_g(x).
\end{align*}
Similarly, for every $t\in\mathbb{R}$,
\begin{align*}
\int_{u^{-1}(t)} f(x)\,d\mathcal{H}^{n-1}_g(x)
=
\sum_{i=1}^{\infty}
\int_{u^{-1}(t)} f_i(x)\,d\mathcal{H}^{n-1}_g(x),
\end{align*}
again by the monotone convergence theorem on the [measure space](/page/Measure%20Space) $(u^{-1}(t),\mathcal{H}^{n-1}_g)$. The chart-supported coarea identity proved below gives the measurability of each function
\begin{align*}
t &\mapsto
\int_{u^{-1}(t)} f_i(x)\,d\mathcal{H}^{n-1}_g(x),
\end{align*}
and the displayed pointwise monotone limit then gives measurability of the corresponding function for $f$. Tonelli's theorem for nonnegative [measurable functions](/page/Measurable%20Functions) then gives
\begin{align*}
\int_{\mathbb{R}}
\left(
\int_{u^{-1}(t)} f(x)\,d\mathcal{H}^{n-1}_g(x)
\right)
\,d\mathcal{L}^1(t)
=
\sum_{i=1}^{\infty}
\int_{\mathbb{R}}
\left(
\int_{u^{-1}(t)} f_i(x)\,d\mathcal{H}^{n-1}_g(x)
\right)
\,d\mathcal{L}^1(t).
\end{align*}
Thus it suffices to prove the formula for a nonnegative Borel function whose support is contained in one coordinate chart.
[/step]
[step:Prove the coarea identity inside one coordinate chart]
Let $(U,\varphi)$ be a coordinate chart, write $V=\varphi(U)\subset\mathbb{R}^n$, and let
\begin{align*}
\psi:V &\to U \\
y &\mapsto \varphi^{-1}(y)
\end{align*}
be the inverse coordinate map. Let $\mathcal{L}^n$ denote $n$-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb{R}^n$, let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$, and let $\mathcal{H}^{n-1}_{\mathbb{R}^n}$ denote $(n-1)$-dimensional [Hausdorff measure](/page/Hausdorff%20Measure) computed with the Euclidean metric on $\mathbb{R}^n$. Assume $h:M\to[0,\infty]$ is Borel measurable and $\operatorname{supp}h\subset U$. Define the coordinate representatives
\begin{align*}
v:V &\to \mathbb{R} \\
y &\mapsto u(\psi(y)),
\end{align*}
and
\begin{align*}
H:V &\to [0,\infty] \\
y &\mapsto h(\psi(y)).
\end{align*}
Since $u$ is Lipschitz with respect to $d_g$ and the Riemannian distance is locally bi-Lipschitz equivalent to the Euclidean distance in a coordinate chart on compact subsets of $V$, the function $v$ is locally Lipschitz. By [Rademacher's Theorem](/page/Rademacher%27s%20Theorem), $v$ is differentiable for $\mathcal{L}^n$-a.e. $y\in V$. Let $D\subset V$ denote the Borel set of differentiability points of $v$. For each $a\in\{1,\dots,n\}$, define $q_a:V\to\mathbb{R}$ to be the Borel representative obtained as the pointwise limit of the coordinate difference quotients on $D$ and $0$ on $V\setminus D$; then $q_a(y)=\partial_{y_a}v(y)$ for every $y\in D$. Define the Borel map
\begin{align*}
q:V &\to \mathbb{R}^n \\
y &\mapsto (q_1(y),\dots,q_n(y)).
\end{align*}
At points $y\in D$, $q(y)=\nabla v(y)$ is the Euclidean gradient of $v$ at $y$.
Let $G(y)=(g_{ab}(y))_{a,b=1}^n$ denote the matrix of the metric tensor in the chart, and define the volume density
\begin{align*}
\omega:V &\to (0,\infty) \\
y &\mapsto \sqrt{\det G(y)}.
\end{align*}
At points where $v$ is differentiable, the Riemannian gradient satisfies
\begin{align*}
|\nabla u(\psi(y))|_g^2
=
\sum_{a,b=1}^n g^{ab}(y)\,\partial_{y_a}v(y)\,\partial_{y_b}v(y),
\end{align*}
where $(g^{ab}(y))_{a,b=1}^n=G(y)^{-1}$. The coordinate formula for Riemannian volume gives
\begin{align*}
d\operatorname{vol}_g(\psi(y))=\omega(y)\,d\mathcal{L}^n(y).
\end{align*}
Define the nonnegative Borel map
\begin{align*}
A:V &\to [0,\infty] \\
y &\mapsto
\begin{cases}
\displaystyle
H(y)\frac{\omega(y)|\nabla u(\psi(y))|_g}{|q(y)|}, & \text{if } y\in D \text{ and } |q(y)|>0,\\
0, & \text{otherwise}.
\end{cases}
\end{align*}
This map is Borel because $H$, $\omega$, $q$, and $y\mapsto |\nabla u(\psi(y))|_g$ on $D$ are Borel after the representative choice above. Because $G(y)$ is positive definite, $|q(y)|=0$ if and only if $|\nabla u(\psi(y))|_g=0$ at every $y\in D$. Hence
\begin{align*}
A(y)|\nabla v(y)|
=
H(y)\omega(y)|\nabla u(\psi(y))|_g
\end{align*}
for $\mathcal{L}^n$-a.e. $y\in V$.
We use the [Euclidean Coarea Formula](/page/Coarea%20Formula) in its open-set locally Lipschitz form. If the cited statement is given only for globally Lipschitz functions on Euclidean space, this open-set form follows by exhausting $V$ by relatively compact open sets, restricting to pieces on which $v$ is Lipschitz, extending those restrictions to globally Lipschitz functions by the McShane [extension theorem](/theorems/59), applying the Euclidean formula there, and passing to the limit by monotone convergence. Its hypotheses are verified: $V\subset\mathbb{R}^n$ is open, $v:V\to\mathbb{R}$ is locally Lipschitz, and $A:V\to[0,\infty]$ is Borel measurable by construction. Applying the formula to $V$, $v$, and $A$ gives
\begin{align*}
\int_V H(y)\omega(y)|\nabla u(\psi(y))|_g\,d\mathcal{L}^n(y)
=
\int_{\mathbb{R}}
\left(
\int_{v^{-1}(t)} A(y)\,d\mathcal{H}^{n-1}_{\mathbb{R}^n}(y)
\right)
\,d\mathcal{L}^1(t).
\end{align*}
We also use the [Area Formula](/page/Area%20Formula) for countably $\mathcal{H}^{n-1}$-rectifiable hypersurfaces. Define the nondifferentiability set and critical set by
\begin{align*}
N&:=V\setminus D, \\
C&:=\{y\in D: |q(y)|=0\}.
\end{align*}
The set $N$ has $\mathcal{L}^n(N)=0$ by [Rademacher's theorem](/page/Rademacher's%20Theorem), and $|\nabla v|=0$ on $C$ in the representative used above. Applying the Euclidean coarea formula to the nonnegative Borel functions $\mathbb{1}_N:V\to[0,\infty]$ and $\mathbb{1}_C:V\to[0,\infty]$ gives
\begin{align*}
\int_{\mathbb{R}}\mathcal{H}^{n-1}_{\mathbb{R}^n}(N\cap v^{-1}(t))\,d\mathcal{L}^1(t)
&=
\int_N |\nabla v(y)|\,d\mathcal{L}^n(y)
=0, \\
\int_{\mathbb{R}}\mathcal{H}^{n-1}_{\mathbb{R}^n}(C\cap v^{-1}(t))\,d\mathcal{L}^1(t)
&=
\int_C |\nabla v(y)|\,d\mathcal{L}^n(y)
=0.
\end{align*}
Hence $\mathcal{H}^{n-1}_{\mathbb{R}^n}((N\cup C)\cap v^{-1}(t))=0$ for $\mathcal{L}^1$-a.e. $t$. On every compact subset of $V$, the smooth positive-definite matrix $G$ has eigenvalues bounded above and below away from $0$, so $\mathcal{H}^{n-1}_g$ and $\mathcal{H}^{n-1}_{\mathbb{R}^n}$ are locally comparable. Exhausting $V$ by compact subsets then gives $\mathcal{H}^{n-1}_g(\psi((N\cup C)\cap v^{-1}(t)))=0$ for the same full-measure set of $t$.
For $\mathcal{L}^1$-a.e. $t\in\mathbb{R}$, the Euclidean coarea theorem also gives that $v^{-1}(t)\setminus (N\cup C)$ is countably $(n-1)$-rectifiable. Fix such a value of $t$ and a point $y\in v^{-1}(t)\setminus (N\cup C)$ at which the approximate tangent hyperplane is $\ker dv_y$. Put
\begin{align*}
\nu(y):=\frac{\nabla v(y)}{|\nabla v(y)|}\in\mathbb{R}^n.
\end{align*}
Choose an Euclidean [orthonormal basis](/page/Orthonormal%20Basis) of $\nu(y)^\perp$ and let $E(y):\mathbb{R}^{n-1}\to\mathbb{R}^n$ be the matrix whose columns are that basis. The squared tangential Jacobian of $\psi$ from the Euclidean hypersurface measure to the Riemannian hypersurface measure is
\begin{align*}
J_{n-1}^{\tau}\psi(y)^2
&=\det\left(E(y)^\top G(y)E(y)\right).
\end{align*}
The block determinant identity for the orthogonal matrix with columns $E(y)$ and $\nu(y)$ gives
\begin{align*}
\det\left(E(y)^\top G(y)E(y)\right)
&=\det G(y)\,\nu(y)^\top G(y)^{-1}\nu(y) \\
&=\omega(y)^2\frac{\nabla v(y)^\top G(y)^{-1}\nabla v(y)}{|\nabla v(y)|^2} \\
&=\left(\frac{\omega(y)|\nabla u(\psi(y))|_g}{|\nabla v(y)|}\right)^2.
\end{align*}
Taking the positive square root gives
\begin{align*}
J_{n-1}^{\tau}\psi(y)
=
\frac{\omega(y)|\nabla u(\psi(y))|_g}{|\nabla v(y)|}.
\end{align*}
For such $t$, the [area formula](/theorems/3075) applied to the smooth diffeomorphism $\psi:v^{-1}(t)\setminus (N\cup C)\to (u^{-1}(t)\cap U)\setminus\psi(N\cup C)$ gives
\begin{align*}
d\mathcal{H}^{n-1}_g(\psi(y))
=
\frac{\omega(y)|\nabla u(\psi(y))|_g}{|\nabla v(y)|}
\,d\mathcal{H}^{n-1}_{\mathbb{R}^n}(y)
\end{align*}
on $v^{-1}(t)\setminus (N\cup C)$. Since $A=0$ on $N\cup C$ and these slices are null for these values of $t$, this implies
\begin{align*}
\int_{v^{-1}(t)} A(y)\,d\mathcal{H}^{n-1}_{\mathbb{R}^n}(y)
=
\int_{u^{-1}(t)\cap U} h(x)\,d\mathcal{H}^{n-1}_g(x)
\end{align*}
for $\mathcal{L}^1$-a.e. $t\in\mathbb{R}$. Substituting this into the preceding identity and using the coordinate formula for volume yields
\begin{align*}
\int_M h(x)|\nabla u(x)|_g\,d\operatorname{vol}_g(x)
=
\int_{\mathbb{R}}
\left(
\int_{u^{-1}(t)} h(x)\,d\mathcal{H}^{n-1}_g(x)
\right)
\,d\mathcal{L}^1(t).
\end{align*}
[guided]
The purpose of this step is to check that the Riemannian formula is exactly the Euclidean formula written in curved coordinates. Let $(U,\varphi)$ be a coordinate chart, let $V=\varphi(U)$, and write
\begin{align*}
\psi:V &\to U \\
y &\mapsto \varphi^{-1}(y).
\end{align*}
Let $\mathcal{L}^n$ denote $n$-dimensional Lebesgue measure on $\mathbb{R}^n$, let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$, and let $\mathcal{H}^{n-1}_{\mathbb{R}^n}$ denote $(n-1)$-dimensional Hausdorff measure computed with the Euclidean metric on $\mathbb{R}^n$. Assume that $h:M\to[0,\infty]$ is Borel measurable and supported in $U$. We transfer all objects to $V$ by defining
\begin{align*}
v:V &\to \mathbb{R} \\
y &\mapsto u(\psi(y)),
\end{align*}
and
\begin{align*}
H:V &\to [0,\infty] \\
y &\mapsto h(\psi(y)).
\end{align*}
The function $v$ is locally Lipschitz because $u$ is Lipschitz on $M$ and, on compact subsets of a coordinate chart, the Riemannian distance and Euclidean coordinate distance are bi-Lipschitz equivalent. [Rademacher's Theorem](/page/Rademacher%27s%20Theorem) therefore gives differentiability of $v$ for $\mathcal{L}^n$-a.e. $y\in V$. Let $D\subset V$ be the Borel set of differentiability points. For each coordinate direction $a\in\{1,\dots,n\}$, the difference quotients defining $\partial_{y_a}v$ are Borel functions wherever they are defined, so their pointwise limit on $D$ has a Borel representative. Denote this representative by $q_a:V\to\mathbb{R}$, set it equal to $0$ on $V\setminus D$, and define
\begin{align*}
q:V &\to \mathbb{R}^n \\
y &\mapsto (q_1(y),\dots,q_n(y)).
\end{align*}
Then $q(y)=\nabla v(y)$ at every $y\in D$. This representative is the reason the coarea weight below is genuinely Borel, rather than merely defined almost everywhere.
Let $G(y)=(g_{ab}(y))_{a,b=1}^n$ be the coordinate matrix of the metric and let
\begin{align*}
\omega:V &\to (0,\infty) \\
y &\mapsto \sqrt{\det G(y)}
\end{align*}
be the Riemannian volume density. The coordinate formula for volume says
\begin{align*}
d\operatorname{vol}_g(\psi(y))=\omega(y)\,d\mathcal{L}^n(y).
\end{align*}
At differentiability points of $v$, if $(g^{ab}(y))_{a,b=1}^n=G(y)^{-1}$, then
\begin{align*}
|\nabla u(\psi(y))|_g^2
=
\sum_{a,b=1}^n g^{ab}(y)\,\partial_{y_a}v(y)\,\partial_{y_b}v(y).
\end{align*}
Now we choose the Euclidean coarea weight so that the Euclidean Jacobian $|\nabla v|$ cancels and leaves the Riemannian integrand. Define the map
\begin{align*}
A:V &\to [0,\infty] \\
y &\mapsto
\begin{cases}
\displaystyle
H(y)\frac{\omega(y)|\nabla u(\psi(y))|_g}{|q(y)|}, & \text{if } y\in D \text{ and } |q(y)|>0,\\
0, & \text{otherwise}.
\end{cases}
\end{align*}
This definition is legitimate because $H$, $\omega$, $q$, and the coordinate formula for $|\nabla u|_g$ are Borel on $D$, while the value on $V\setminus D$ is fixed to be $0$. Since the metric matrix $G(y)$ is positive definite, $q(y)=0$ exactly when $\nabla u(\psi(y))=0$ at differentiability points. Thus
\begin{align*}
A(y)|\nabla v(y)|
=
H(y)\omega(y)|\nabla u(\psi(y))|_g
\end{align*}
holds for $\mathcal{L}^n$-a.e. $y\in V$.
We now apply the [Euclidean Coarea Formula](/page/Coarea%20Formula) in its open-set locally Lipschitz form. If the referenced theorem is stated only for globally Lipschitz maps on Euclidean space, we obtain this form by exhausting $V$ with relatively compact open pieces, using local Lipschitzness of $v$ on each piece, extending each restricted function to a globally Lipschitz function by the McShane extension theorem, and passing through monotone convergence. Its hypotheses are verified: $V\subset\mathbb{R}^n$ is open, $v:V\to\mathbb{R}$ is locally Lipschitz, and $A:V\to[0,\infty]$ is Borel measurable by construction. Therefore
\begin{align*}
\int_V A(y)|\nabla v(y)|\,d\mathcal{L}^n(y)
=
\int_{\mathbb{R}}
\left(
\int_{v^{-1}(t)} A(y)\,d\mathcal{H}^{n-1}_{\mathbb{R}^n}(y)
\right)
\,d\mathcal{L}^1(t).
\end{align*}
Using the preceding identity for $A|\nabla v|$, the left-hand side becomes
\begin{align*}
\int_V H(y)\omega(y)|\nabla u(\psi(y))|_g\,d\mathcal{L}^n(y).
\end{align*}
It remains to identify the level-set measure. This is not a purely formal coordinate rewrite: for a Lipschitz function, level sets are hypersurfaces only in the rectifiable a.e. sense, and critical level sets may be large at exceptional values of $t$. Define the nondifferentiability and critical sets by
\begin{align*}
N&:=V\setminus D, \\
C&:=\{y\in D: |q(y)|=0\}.
\end{align*}
We must discard both sets. The critical set $C$ is where the coarea weight was defined to vanish because the gradient is zero; the nondifferentiability set $N$ is where the coordinate gradient was not defined and where we also set $A$ equal to $0$. [Rademacher's theorem](/theorems/3069) gives $\mathcal{L}^n(N)=0$, and $|\nabla v|=0$ on $C$ in the chosen representative. Applying the same Euclidean coarea formula to the nonnegative Borel maps $\mathbb{1}_N:V\to[0,\infty]$ and $\mathbb{1}_C:V\to[0,\infty]$ gives
\begin{align*}
\int_{\mathbb{R}}\mathcal{H}^{n-1}_{\mathbb{R}^n}(N\cap v^{-1}(t))\,d\mathcal{L}^1(t)
&=
\int_N |\nabla v(y)|\,d\mathcal{L}^n(y)
=0, \\
\int_{\mathbb{R}}\mathcal{H}^{n-1}_{\mathbb{R}^n}(C\cap v^{-1}(t))\,d\mathcal{L}^1(t)
&=
\int_C |\nabla v(y)|\,d\mathcal{L}^n(y)
=0.
\end{align*}
Because both integrands on the left are nonnegative, $\mathcal{H}^{n-1}_{\mathbb{R}^n}((N\cup C)\cap v^{-1}(t))=0$ for $\mathcal{L}^1$-a.e. $t$. This does not say that every bad level is null; it says exactly what is needed for an identity integrated over $t$. On compact subsets of $V$, the eigenvalues of the smooth positive-definite matrix $G$ are bounded above and below away from $0$, so $\mathcal{H}^{n-1}_g$ and $\mathcal{H}^{n-1}_{\mathbb{R}^n}$ are locally comparable. Exhausting $V$ by compact subsets transfers the zero-measure conclusion to $\mathcal{H}^{n-1}_g$ on $\psi((N\cup C)\cap v^{-1}(t))$ for the same full-measure set of levels.
For those values of $t$, the Euclidean coarea theorem also gives that $v^{-1}(t)\setminus (N\cup C)$ is countably $(n-1)$-rectifiable. The coordinate change $\psi$ sends $v^{-1}(t)\setminus (N\cup C)$ to $(u^{-1}(t)\cap U)\setminus\psi(N\cup C)$. We now compute the tangential Jacobian. At an approximate tangent point $y\in v^{-1}(t)\setminus (N\cup C)$, the Euclidean tangent hyperplane is $\ker dv_y$, whose unit normal is
\begin{align*}
\nu(y):=\frac{\nabla v(y)}{|\nabla v(y)|}.
\end{align*}
Choose an Euclidean orthonormal basis of $\nu(y)^\perp$, and let $E(y):\mathbb{R}^{n-1}\to\mathbb{R}^n$ be the matrix with those basis vectors as columns. The pullback of the Riemannian metric to this tangent plane has Gram matrix $E(y)^\top G(y)E(y)$, so
\begin{align*}
J_{n-1}^{\tau}\psi(y)^2
&=\det\left(E(y)^\top G(y)E(y)\right).
\end{align*}
If we append $\nu(y)$ as the last column, the resulting matrix is orthogonal. The Schur-complement block determinant identity gives
\begin{align*}
\det\left(E(y)^\top G(y)E(y)\right)
&=\det G(y)\,\nu(y)^\top G(y)^{-1}\nu(y) \\
&=\omega(y)^2\frac{\nabla v(y)^\top G(y)^{-1}\nabla v(y)}{|\nabla v(y)|^2} \\
&=\left(\frac{\omega(y)|\nabla u(\psi(y))|_g}{|\nabla v(y)|}\right)^2.
\end{align*}
The Jacobian is nonnegative, so taking square roots gives
\begin{align*}
J_{n-1}^{\tau}\psi(y)
=
\frac{\omega(y)|\nabla u(\psi(y))|_g}{|\nabla v(y)|}.
\end{align*}
The [Area Formula](/page/Area%20Formula) therefore gives
\begin{align*}
d\mathcal{H}^{n-1}_g(\psi(y))
=
\frac{\omega(y)|\nabla u(\psi(y))|_g}{|\nabla v(y)|}
\,d\mathcal{H}^{n-1}_{\mathbb{R}^n}(y)
\end{align*}
on $v^{-1}(t)\setminus (N\cup C)$. Since $A=0$ on $N\cup C$, and since these discarded slices have zero measure for these values of $t$ with respect to both level measures, we may ignore exactly those slices in the inner integrals. Hence
\begin{align*}
\int_{v^{-1}(t)} A(y)\,d\mathcal{H}^{n-1}_{\mathbb{R}^n}(y)
=
\int_{u^{-1}(t)\cap U} h(x)\,d\mathcal{H}^{n-1}_g(x)
\end{align*}
for $\mathcal{L}^1$-a.e. $t$. Substitution gives the chart-supported identity
\begin{align*}
\int_M h(x)|\nabla u(x)|_g\,d\operatorname{vol}_g(x)
=
\int_{\mathbb{R}}
\left(
\int_{u^{-1}(t)} h(x)\,d\mathcal{H}^{n-1}_g(x)
\right)
\,d\mathcal{L}^1(t).
\end{align*}
[/guided]
[/step]
[step:Sum the chart identities to obtain the global formula]
Apply the chart-supported identity from the previous step to each $f_i=\rho_i f$. Since $\operatorname{supp}f_i\subset U_i$, we have
\begin{align*}
\int_M f_i(x)|\nabla u(x)|_g\,d\operatorname{vol}_g(x)
=
\int_{\mathbb{R}}
\left(
\int_{u^{-1}(t)} f_i(x)\,d\mathcal{H}^{n-1}_g(x)
\right)
\,d\mathcal{L}^1(t)
\end{align*}
for every $i\in\mathbb{N}$. Summing over $i$ and using the monotone-convergence decompositions established in the first step gives
\begin{align*}
\int_M f(x)|\nabla u(x)|_g\,d\operatorname{vol}_g(x)
=
\int_{\mathbb{R}}
\left(
\int_{u^{-1}(t)} f(x)\,d\mathcal{H}^{n-1}_g(x)
\right)
\,d\mathcal{L}^1(t).
\end{align*}
This proves the formula for every nonnegative Borel measurable $f:M\to[0,\infty]$.
[/step]
[step:Pass from Borel functions to completed measurable functions]
Let $f:M\to[0,\infty]$ be measurable with respect to the completion of the Borel $\sigma$-algebra for $\operatorname{vol}_g$. Choose a nonnegative Borel function $F:M\to[0,\infty]$ such that $F=f$ $\operatorname{vol}_g$-a.e., and choose a Borel set $E\subset M$ with $\operatorname{vol}_g(E)=0$ containing the exceptional set $\{x\in M:F(x)\neq f(x)\}$. The left-hand sides agree because $f|\nabla u|_g=F|\nabla u|_g$ $\operatorname{vol}_g$-a.e. Applying the already-proved Borel coarea identity to $\mathbb{1}_E:M\to[0,\infty]$ gives
\begin{align*}
\int_{\mathbb{R}}
\mathcal{H}^{n-1}_g(E\cap u^{-1}(t))
\,d\mathcal{L}^1(t)
=
\int_M \mathbb{1}_E(x)|\nabla u(x)|_g\,d\operatorname{vol}_g(x)
=0.
\end{align*}
Since the integrand is nonnegative, $\mathcal{H}^{n-1}_g(E\cap u^{-1}(t))=0$ for $\mathcal{L}^1$-a.e. $t$. For those values of $t$, the inner level-set integrals of $f$ and $F$ agree. Therefore the coarea identity for $F$ implies the same identity for $f$.
[/step]
[step:Extend the identity from nonnegative functions to integrable signed functions]
Let $f:M\to\mathbb{R}$ be measurable with respect to the completed Riemannian volume measure and suppose the corresponding coarea identity for $|f|$ is finite, meaning that either, and therefore both by the nonnegative formula, of the following quantities is finite:
\begin{align*}
\int_M |f(x)|\,|\nabla u(x)|_g\,d\operatorname{vol}_g(x),
\qquad
\int_{\mathbb{R}}
\left(
\int_{u^{-1}(t)} |f(x)|\,d\mathcal{H}^{n-1}_g(x)
\right)
\,d\mathcal{L}^1(t).
\end{align*} Define the positive and negative parts
\begin{align*}
f^+:M &\to [0,\infty] \\
x &\mapsto \max\{f(x),0\},
\end{align*}
and
\begin{align*}
f^-:M &\to [0,\infty] \\
x &\mapsto \max\{-f(x),0\}.
\end{align*}
Then $f=f^+-f^-$ and $|f|=f^++f^-$. Applying the nonnegative coarea formula to $f^+$ and $f^-$ gives two finite identities; finiteness follows from applying the already-proved nonnegative formula to $|f|=f^++f^-$. Subtracting the identity for $f^-$ from the identity for $f^+$ is legitimate because finiteness of the formula for $|f|$ excludes the indeterminate expression $\infty-\infty$. Therefore
\begin{align*}
\int_M f(x)|\nabla u(x)|_g\,d\operatorname{vol}_g(x)
=
\int_{\mathbb{R}}
\left(
\int_{u^{-1}(t)} f(x)\,d\mathcal{H}^{n-1}_g(x)
\right)
\,d\mathcal{L}^1(t).
\end{align*}
This is precisely the asserted extension to integrable signed functions.
[/step]