I understand the reason why the gradient vector is always orthogonal to the level sets of f, but I just cannot find any notes saying the gradient vector is tangent to the surface.
But it seems to be reasonable if I imagine that when climbing the hill, the steepest path is actually tangent to the hill!!!
So I am now a bit confusing...
On
Let $f:\mathbb R^n \rightarrow \mathbb R$ be the function in question.
I assume that by "the surface", you mean the (hyper)surface in $\mathbb R^{n+1}$ that's defined $$S=\left \{ (x, y)\in\mathbb R^{n+1} \text { such that } y=f(x)\right\}$$
If so, then the normal to S at the point $(x, y) \in \mathbb R^{n+1}$ is given by $$(\nabla f(x), -1)$$
To see this, define $g:\mathbb R^{n+1} \rightarrow \mathbb R$ to be equal to $g(x, y)=y-f(x)$. Then $S$ is the level set of $g$ defined by $g(x,y)=0$. As you know, the gradient of $g$ gives the normal to its level sets, so the normal to $S$ is given by $(-\nabla f(x), 1)$. Multiply that vector by $-1$ and it's still a normal vector.
On
“The hill” is the surface $z=f(x,y)$ in $xyz$-space, but what the gradient points out is the direction of steepest ascent on “the map”, namely the $xy$-plane.
Remember that the gradient of a two-variable function is a vector with only two components, so it doesn't even make sense to talk about how it's pointing in $\mathbf{R}^3$.
You cannot any notes saying the gradient vector is tangent to the surface, because it is not. If you define a surface $S$ as a level curve:$$S=\bigl\{(x,y,z)\in\mathbb R^3\mid f(x,y,z)=c\bigr\},$$and if $p\in S$, then $\nabla f(p)$ is orthogonal to $S$, not tangent to it, since it gives you the direction in which $f$ is groing fastar and, within $S$, $f$ does not grow at all.