Why is the gradient normal? This is a somewhat long discussion so please bear with me.  There is a theorem that I have always been curious about from an intuitive standpoint and that has been glossed over in most textbooks I have read.  Quoting Wikipedia, the theorem is:

The gradient of a function at a point is perpendicular to the level set of $f$ at that point.

I understand the Wikipedia article's proof, which is the standard way of looking at things, but I see the proof as somewhat magical.  It gives a symbolic reason for why the theorem is true without giving much geometric intuition.
The gradient gives the direction of largest increase so it sort of makes sense that a curve that is perpendicular would be constant.  Alas, this seems to be backwards reasoning.  Having already noticed that the gradient is the direction of greatest increase, we can deduce that going in a direction perpendicular to it would be the slowest increase.  But we can't really reason that this slowest increase is zero nor can we argue that going in a direction perpendicular to a constant direction would give us a direction of greatest increase.
I would also appreciate some connection of this intuition to Lagrange multipliers which is another somewhat magical theorem for me.  I understand it because the algebra works out but what's going on geometrically?
Finally, what does this say intuitively about the generalization where we are looking to: maximize $f(x,y)$ where $g(x,y) > c$.
I have always struggled to find the correct internal model that would encapsulate these ideas.
 A: If gradient is not perpendicular to the level curve, it will have some component along the level curve. This means that function's value will increase if you move in that direction, but on a level curve, function's value can not increase (or decrease) so gradient can not have component along level curve and this is possible only when gradient is perpendicular to level curve.
A: A very useful way to think about the gradient (or more generally, the first derivative of any function on some Euclidean space) is as "the thing that gives you the best linear approximation to the function at a given point." More explicitly, If p is some point in space, then if you want to approximate f by a linear function near p in the best way possible, then you should take your linear function to be:
L(x) = f(p) + f'(p).(x-p), 
where f'(p) is the gradient at p. (That generalizes, by the way, to vector-valued functions: then f'(p) is a linear transformation, not simply a vector.) The point is that if you look at your function very closely near the point p, then it looks more and more like that linear function, and the approximation just gets better as you look closer. In particular, that linear approximation completely captures both the direction of the level curve of f through p and the direction of fastest growth. Now you should visualize a linear function of two variables, whose graph is simply a slanted plane, and it should be obvious that the level curves (which are horizontal lines embedded in the plane) are perpendicular to the slant of the plane, if you're visualizing it correctly. 
One subtlety that my explanation doesn't completely cover is what happens when the gradient is zero: that would give you the linear approximation 
L(x) = f(p),
i. e., constant. Of course, the function f may not literally be constant near p, so you might still want to know what the level set looks like and where the fastest growth is. To get that information, you will have to use the higher derivatives. All the gradient tells you is that:
(1) the fastest growth, wherever it is, is slower than any linear function if you're close enough to p, and 
(2) the level set through p is "probably" experiencing some sort of singularity there, like self-intersection, or possibly degenerating to a single point (in the case of a local extremum).  
A: Start in the 1-variable case. The derivative of a function at a point is a number that can be considered a vector that either points left or right. It is by default perpendicular to the level curve which is a point. That vector, f'(c), is combined to form
(f'(c),-1) which is a vector perpendicular to the graph of the function y=f(x) at the point(c,f(c)). 
For a function z=g(x,y) of two variables and for most points, there is a neighborhood of the point at which the level curve g(c,d)=K can be written as the graph of a function y=f(x) --- that is a local thing. So the perpendicular to the tangent of that graph (f'(c),-1) is up to a constant the gradient of g. 
Finally, consider the case of a linear function. z=Ax+By. The vector (A,B) is perpendicular to the line Constant =Ax +By.
A: The proof I usually see: Choose an arbitrary unit length tangent vector on the level set, and write it with coordinates.  If you take the inner product of this vector with the gradient, the sum you get is the definition of the directional derivative along that vector, and therefore zero.
I think there is a more geometric way to think of it by taking a linear approximation to your function near a noncritical point, and restricting the linear function to a sphere.  After a suitable rotation, the extrema lie at poles, and the level set lies at the equator.
A: This is essentially Andrew Stacey's answer, but a bit lower level. This is the story I actually try to get my calculus 3 students to understand.
Let $F: \mathbb{R}^2 \to \mathbb{R}$.  Then the derivative $D_{F,p}$ is a linear map from $D_{F,p}:\mathbb{R}^2 \to \mathbb{R}$, whose matrix with respect to the standard basis is $[ \dfrac{\partial F}{\partial x} \dfrac{\partial F}{\partial y}]$.
This is the unique linear map which satisfies $F(p+h) = F(p)+D_{F,p}(h)+Error(h)$, where $\displaystyle\lim_{h \to 0} \dfrac{|Error(h)|}{|h|} = 0$.  Notice that $p$ and $h$ are both vectors in $\mathbb{R}^2$.
The cool thing about linear maps from $\mathbb{R}^n \to \mathbb{R}$ is that they look like dot products!  In this case, with $h = \langle a,b \rangle$, evaluating the derivative at point $p$, then we have $D_{F,p}(\langle a,b \rangle) = [ \dfrac{\partial F}{\partial x} \dfrac{\partial F}{\partial y}] \begin{bmatrix}a \\\\ b\end{bmatrix} = \dfrac{\partial F}{\partial x}a + \dfrac{\partial F}{\partial y}b = \langle \dfrac{\partial F}{\partial x} \dfrac{\partial F}{\partial y} \rangle \cdot \langle a, b\rangle$,
(with the partials evaluated at $p$). This alternative viewpoint on the derivative is useful, because it gives a different geometric interpretation of the derivative.  We call $\langle \dfrac{\partial F}{\partial x} \dfrac{\partial F}{\partial y} \rangle$ the gradient of $F$.
Now we are interested in the curve $F(x,y) = 0$.  Given a point $p=(x_1,y_1)$ on this curve, the tangent direction will be the vector $h$ for which $D_{F,p}(h) = 0$, because to stay on the curve, the value of the function should not change to first order. Using the geometric interpretation in terms of dot products, we can see that $\langle \dfrac{\partial F}{\partial x} \dfrac{\partial F}{\partial y} \rangle \cdot \langle h_1, h_2\rangle = 0$, or geometrically that the gradient is perpendicular to the tangent direction!
A: If you are standing on a level set and want to walk some small distance d and get as far as possible from the level set, you want to walk along the normal. Otherwise, if the path you take has a tangent component, it will tend to keep you closer to the level set if d is small enough compared to the size of the level set. Furthermore, getting as far as possible from your level set is approximately the same as walking to the highest/lowest level curve in range, with the approximation improving as d shrinks.
A: Imagine two surfaces: $S_1 : f(x,y,z) = c_1$ and $S_2 : f(x,y,z) = c_2$ where $c_1 < c_2 $. Now let's see what happens in a region $ x \in (a,a+\Delta x),y \in (b,b+\Delta y) $ in which $S_1$ is approximately flat. Let's take the perpendicular to $S_1$ lines that extend from the corner points of $S_1$ to $S_2$. In the region that these lines define, $S_2$ is going to have a point $P_{min}$ with minimal distance from $S_1$ and a point $P_{max}$ with maximal distance from $S_1$. Let's define as $D$ the difference between their respective distances from $S_1$. Now we can arbitrarily choose a $D' < D$ and a surface of the same form will exist between $S_1$ and $S_2$ with a point $P'_{max}$ with maximal distance $D'$ from $S_1$. This is allowed because $f(x,y,z)$ is continuous and the surfaces cannot intersect. Let $P'_{min}$ be the point with minimal distance from $S_1$. Then the difference of their respective 
distances will be less than $D'$. If we choose a $D'$ small enough the resulting surface $S' : f(x,y,z) = c'$ will be approximately flat and 
parallel to $S_1$ in the specified region. Now the gradient $\nabla f$ points to the direction of the biggest increase of $f(x,y,z)$. Inside the region, the shortest path from any point on $S_1$ to $S'$ is through the perpendicular line. Since the same will happen for all surfaces between $S_1$ and $S'$, the gradient will be perpendicular to the surface.
A: $\newcommand{\+}{^{\dagger}}%
 \newcommand{\angles}[1]{\left\langle #1 \right\rangle}%
 \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}%
 \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}%
 \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}%
 \newcommand{\dd}{{\rm d}}%
 \newcommand{\down}{\downarrow}%
 \newcommand{\ds}[1]{\displaystyle{#1}}%
 \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}%
 \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}%
 \newcommand{\fermi}{\,{\rm f}}%
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}%
 \newcommand{\half}{{1 \over 2}}%
 \newcommand{\ic}{{\rm i}}%
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}%
 \newcommand{\isdiv}{\,\left.\right\vert\,}%
 \newcommand{\ket}[1]{\left\vert #1\right\rangle}%
 \newcommand{\ol}[1]{\overline{#1}}%
 \newcommand{\pars}[1]{\left( #1 \right)}%
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}%
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}%
 \newcommand{\sech}{\,{\rm sech}}%
 \newcommand{\sgn}{\,{\rm sgn}}%
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}%
 \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$
Take two points $\vec{R}$ and $\vec{R} + \delta\vec{R}$ in a surface defined by $\fermi\pars{\vec{r}} = 0$. Then, $\fermi\pars{\vec{R}} = 0$ and $\fermi\pars{\vec{R} + \delta\vec{R}} = 0$. Then,
$$
0
= \lim_{\verts{\delta\vec{R}} \to 0}
{\fermi\pars{\vec{R} + \delta\vec{R}} -\fermi\pars{\vec{R}}
  \over \verts{\delta\vec{R}}} 
=
\nabla_{\vec{R}}\fermi\pars{\vec{R}}\cdot\lim_{\verts{\delta\vec{R}} \to 0}{\delta\vec{R} \over \verts{\delta\vec{R}}}
$$ 
Provided that
$\hat{n}\pars{\vec{R}}
 \equiv \lim_{\verts{\delta\vec{R}} \to 0}\pars{%
{\delta\vec{R}/ \verts{\delta\vec{R}}}}$ exists, $\hat{n}\pars{\vec{R}}$ is a 'tangent vector' to the surface at point $\vec{R}$ which yields
$\ds{\nabla_{\vec{R}}\fermi\pars{\vec{R}}\cdot\hat{n}\pars{\vec{R}} = 0}$. Then
$\ds{\color{#00f}{\large\nabla_{\vec{R}}\fermi\pars{\vec{R}}\perp\hat{n}\pars{\vec{R}}}}$. 
A: As for Lagrange multipliers (taking the geometric interpretation of the gradient as given), I think the best way to see what's going on is to consider the 2-variable case.  The idea is that you're generally constrained to some curve g(x,y)=0 (which is just a particular level curve of the function g(x,y)).  Consider the level curves of the objective function f(x,y).  If, moving along the constraint curve, we happen to be crossing one of these level curves of f(x,y), then (subject to the constraint) we are either increasing or decreasing the value of f(x,y).  Therefore, we cannot be at a local extremum.  So, a local extremum must occur where g(x,y)=0 is tangent to a level curve of f(x,y).  This is equivalent to the condition that vectors perpendicular to the curves point in the same direction: grad(f)=λgrad(g).
A: Kim, remember the level sets' names in practical applications are isoclines, isotherms, isobars, isoutility curves: they are loci of places which all share the same (scalar) value, say 3000 feet elevation.
 (source)
$$\text{ level set of 3000 feet } \overset{def}{=} \ \{\ \forall \text{ points } p \text{ in the surrounding area }  \ | \ \mathtt{height}(p) = 3000 \ \}$$
(In the picture above the level curves are continuous)
The shortest path between $p$ on $\text{the level set of 3000 feet}$ and (somewhere on) $\text{the level set of 4000 feet}$ will also be the fastest way or the least circuitous way up.

If you take this concept to the infinitesimal limit of slightly-higher level sets you get the gradient at $p$.
Alternately you can imagine that someone cut a walkable path on the side of the mountain. It's walkable because it's completely flat, i.e. it follows a level curve (constant height). Walking north near my red $p$ look to your left. You are looking straight uphill. This is the gradient direction.

Now keep in mind this is just one gradient at one point. Maybe you heard of a gradient field which sticks an uphill arrow at every point in the surrounding area. If you can envisage the gradient field of the pictures above (from volcano.oregonstate.edu) then you'll know you've got the concept.
edit: The original picture I was linking to moved; here are some more topographic maps.

(source)

(source: Wayback Machine)

(source)
(source: Wayback Machine)

(source)

(source)
A: The gradient of a function is normal to the level sets because it is defined that way.  The gradient of a function is not the natural derivative.  When you have a function $f$, defined on some Euclidean space (more generally, a Riemannian manifold) then its derivative at a point, say $x$, is a function $d_xf(v)$ on tangent vectors.  The intuitive way to think of it is that $d_xf(v)$ answers the question:

If I move infinitesimally in the direction $v$, what happens to $f$?

So $d_xf(v)$ is not itself a tangent vector.  However, as we have an inner product lying around, we can convert it into a tangent vector which we call $\nabla f$.  This represents the question:

What tangent vector $u$ at $x$ best represents $d_xf(v)$?

What we mean by "best represents" is that $u$ should satisfy the condition:

$\langle u,v\rangle  = d_xf(v)$ for all tangent vectors $v$.

Now we look at the level set of $f$ through $x$.  If $v$ is a tangent vector at $x$ which is tangent to the level set then $d_xf(v) = 0$ since $f$ doesn't change if we go (infinitesimally) in the direction of $v$.  Hence our vector $\nabla f$ (aka $u$ in the question) must satisfy $\langle\nabla f, v\rangle = 0$.  That is, $\nabla f$ is normal to the set of tangent vectors at $x$ which are tangent to the level set.
For a generic $x$ and a generic $f$ (i.e. most of the time), the set of tangent vectors at $x$ which are tangent to the level set of $f$ at $x$ is codimension $1$ so this specifies $\nabla f$ up to a scalar multiple.  The scalar multiple can be found by looking at a tangent vector $v$ such that $f$ does change in the $v$-direction.  If no such $v$ exists, then $\nabla f = 0$, of course.
A: I like this intuition:
The bundle of tangent vectors to the surface at a point live in the tangent plane at that point.  The tangent[s] to the level set at that point are exactly the vectors in the tangent plane whose "vertical" component is zero.  The vector[s] pointing in the direction of greatest increase are those with the largest relative "vertical" components.  Plane geometry (in the tangent plane) shows that these must to be perpendicular.
