This is essentially asking again the following question: Formula for the induced Hessian, which never received an answer. Below I record my thoughts.
For $\Sigma^{n-1}\subset M^n$, I want to show $$|\nabla^2_\Sigma u|^2=|\nabla^2 u|^2-2|\nabla|\nabla u||^2+|\nabla^2 u(\nu,\nu)|^2.$$ Since this is an equality of functions, it suffices to show this equality pointwise. Therefore, we use geodesic normal coordinates, with $e_n=\nu$. The LHS in coordinates becomes $$\text{LHS}=\sum_{i,j}^{n-1}u_{ij}^2,$$ while the RHS in coordinates (I hope) becomes \begin{align*} \text{RHS}&=\left(\sum_{i,j}^nu_{ij}^2\right)-2\sum_{i=1}^nu_{in}+u_{nn}^2\\ &=\left(\sum_{i,j=1}^{n-1}u_{ij}^2+2\sum_{i=1}^{n-1}u_{in}^2+u_{nn}^2\right)-2\sum_{i=1}^nu_{in}+u_{nn}^2\\ &=\text{LHS} \end{align*} Thus, the only term requiring justification is why in geodesic normal coordinates, $$|\nabla|\nabla u||^2=\sum_{i=1}^n u_{in}^2.$$ We compute this term by using $$|\nabla v|^2= \sum_{i=1}^n v_i^2,$$ with $v=|\nabla u|$, and also $$|\nabla u|_i=\sum_{j=1}^n \frac{u_{j}u_{ij}}{|\nabla u|}.$$ Altogether, we must show $$\sum_{i=1}^n\left(\sum_{j=1}^n \frac{u_ju_{ij}}{|\nabla u|}\right)^2=\sum_{i=1}^n u_{in}^2.$$ Is there a mistake somewhere, or can someone help me finish the last step in the computation?
Got it! The answer relies on a simple calc 3 observation.
The computation finishes upon realizing that for the level set is perpendicular to the gradient, so $|\nabla u|=u_n$ (i.e. $u_j=0$ for $j=1,...,n-1$). Therefore,
\begin{align*} \sum_{i=1}^n\left(\sum_{j=1}^n \frac{u_ju_{ij}}{|\nabla u|}\right)^2&=\sum_{i=1}^n\left(\frac{u_nu_{in}}{|\nabla u|}\right)^2\\ &=\sum_{i=1}^nu_{in}^2. \end{align*}
If anyone can produce a coordinate-free proof, I'd be happy to accept your answer. In particular, I would be interested to see if there's a way to express that $\Sigma$ is perpendicular to $\nabla u$ without coordinates.