Let $M=f^{-1}(\{0\})$ be a $d$ dimensional manifold situated in $\mathbb{R}^D$, where $D>d$. We assume $f:\mathbb{R}^D\to \mathbb{R}^p$ where $p,D,d$ are related by $p=D-d$. For those interested in regularity, we assume $f$ is $C^\infty$ and its Jacobian has full rank everywhere on $M$.
We are trying to find explicit formulas for the orthogonal projection matrix from $\mathbb{R}^D$ onto the tangent space of $M$ at a given point $x\in M$. Explicitly we have the following
Question: What are explicit formulas for the orthogonal projection matrix $P: \mathbb{R}^D \to \mathbb{R}^{D\times D}$ where $P(x)P(x)=P(x)^2=P(x)$.
What I know: If $d=D-1$, i.e. $p=1$, then the following form is rather well known: $$P(x)=I-n(x)n(x)^T,$$ where $n$ is the normal vector to $M$, $$n(x)=\frac{\nabla f(x)}{\|\nabla f(x)\|}.$$
One way to derive this is is to write each point of $M$ as $x=(y, F(y))$ where $y\in \mathbb{R}^{D-1}$ and $F: \mathbb{R}^{D-1}\to \mathbb{R}$. The relationship between $f$ and $F$ is as follows: $$f^i(x_1,\dotsc, x_D)=F^i(x_1,\dotsc, x_d)-x^{d+i}.$$ Then the intrinsic metric tensor is $g(y)=I+J_F^T J_F$ where $I$ is $d\times d$ identity matrix and $J_F$ is the $1\times d$ Jacobian matrix of $F$. Then, $$\tilde{P}(y) = (I, J_F^T)^T g^{-1}(y) (I, J_F^T),$$ defines an orthogonal projection via $P(x)=\tilde{P}(\pi(x))$ where $\pi$ is the projection of the first $d$ coordinates. Writing out the computations explicitly yields the above form of $P$, it is only a little tedious.
Unfortunately, I am having trouble generalizing this for $d=D-p$, $p>1$. My guess is that we should have something analogous like $P=I-J_f^T J_f$ (note $J_f$ is $p\times D$ so this product, with $I$ a $D\times D$ identity, at least gives the right dimensions of $D\times D$). But I've checked for some examples, it is not correct. We are clearly missing the normalization factor present in the $I-nn^T$ hypersurface case. I apologize that my vector calculus/differential geometry is a bit weak, so if this is well known, I'd gladly take references to read into.
Update: Thought about it a little and realized $P=I_{D\times D}-A^TA$, where $A$ is $p\times D$ is an orthogonal projection if and only if $AA^T=I_{p\times p}$. Since $A$ is rectangular, this means that the rows of $A$ must be orthonormal. So for $A=J_f$ this means we need to orthonormalize the rows. Then we should have our expression for $P$. I'm going to try to work out all the details. Comments/corrections on this approach are welcome.