I would like to get some intuitive feeling for the mean curvature. The mean curvature of a hypersurface in a Riemannian manifold by definition is the trace of the second fundamental form.
Is there any good picture to have in mind, when dealing with it?
Do this experiment: draw a curve in 2D and about some point on that curve, draw a unit normal vector field. Now convince yourself that if you push the curve out along the normal field by a distance $\epsilon$, the length changes by a factor of $1+\epsilon\kappa$ where $\kappa$ is the curvature.
Now consider a 2D surface in 3D. Choose local coordinates so that the first two basis vectors are tangent to the surface and the curvatures along the first two coordinates are the principal curvatures. The infinitesimal change in surface area is $(1+\epsilon\kappa_1)(1+\epsilon\kappa_2)$ -- actually this product is the Jacobian of the normal map, $N_\epsilon$ that pushes the surface out along the normal field a distance epsilon.
In general, for co-dimension one surfaces, the Jacobain of the normal map $N_\epsilon$ is $\Pi_{i=1}^{n-1} (1+\epsilon\kappa_i)$. We can integrate this over the original surface to get the new, n-1 volume of the pushed surface. That is:
$\mathcal{H}^{n-1}(N_\epsilon(W)) = \int_{W} \Pi_{i=1}^{n-1} (1 + \epsilon\kappa_i) d\mathcal{H}^{n-1}$
$ \hspace{1in} = \int_{W} 1 d\mathcal{H}^{n-1} + \epsilon \int_{W} \sum_i \kappa_i d\mathcal{H}^{n-1} + ... + \epsilon^k \int_{W} \sum_{s\in S(k)}\Pi_{i\in s}\kappa_i d\mathcal{H}^{n-1}$
$ \hspace{1.2in} + ... + \epsilon^{n-1} \int_{\partial W} \Pi_{i=1}^{n-1} \kappa_i d\mathcal{H}^{n-1}$
Now we note that the first order term is the integral of the mean curvature. That is, to first order, the change in surface volume is given by the mean curvature. Note that I am using the term mean curvature for what is sometimes called total mean curvature, $\sum_{i=1}^{n-1} \kappa_i$.
(This works in co-dimension $k>1$, but then, because the set of normal directions at any point is $k$-dimensional, it is a bit more involved to get a result that looks like the result here.)
First, rather than viewing the mean curvature as "the trace of the second fundamental form," it might be more intuitive to view it as the average of the principle curvatures, i.e., $$H = \frac{1}{n} \sum_{i=1}^n \kappa_i \;.$$
Second, it may be that the following view of minimal surfaces (those of mean curvature zero) from D. Hoffman and W. H. Meeks III, in their paper, "Minimal surfaces based on the catenoid" [Amer. Math. Monthly 97(8) (1990), 702-730] (ACM link), might help:
“Loosely speaking, one imagines the surface as made up of very many rubber bands, stretched out in all directions; on a minimal surface the forces due to the rubber bands balance out, and the surface does not need to move to reduce tension.”
You can see this accords with Peter Michor's more abstract formulation.
It gives you a normal field along the hyper surface which describes how you should move the surface to shrink the induced surface volume form. This works also for higher codimension.
EDIT: An explicit formula which generalized the answer of Otis Chodosh is the following:
Let $(N,\bar g)$ be a Riemannian manifold and let $M$ be a compact manifold with $\dim(M)<\dim(N)$. Let $f(t,\quad)\in Imm(M,N)$ be a smooth curve of immersions, then $f_t = \partial_t|_0 f\in \Gamma(f^*TN)$ splits into parts $f_t=Tf.f_t^\top + f_t^\bot$ tangent and normal to $M$. Let also '$g=f^\star\bar g$ be the induced metric. Then the variation of the mapping $$ Imm(M,N) \to \Gamma(vol(M)),\qquad f \mapsto vol(g)=vol(f^*\bar g) $$ is given by $$ \partial_t|_0\, vol\Big(f(t,\quad)^\star\bar g\Big) = Tr^g\big(\bar g(\nabla f_t,Tf)\big) vol(g)= \Big({div}^{g}(f_t^{\top})-\bar g\big(f_t^{\bot},Tr^g(S)\big)\Big) vol(g). $$ Here $S$ is the second fundamental form $S\in\Gamma(L^2_{sym}(TM;TM^\bot)$. So if the variation field $f_t$ is always normal to $f(M)$, the first summand vanishes.
For a proof see 5.7 in: Martin Bauer, Philipp Harms, Peter W. Michor: Sobolev metrics on shape space of surfaces. Journal of Geometric Mechanics 3, 4 (2011), 389-438. (pdf)
Perhaps a more geometric way of viewing mean curvature (at least for a hypersurface) is as follows:
If $\varphi: \Sigma^n\hookrightarrow (M^{n+1},g)$ is an embedded (oriented) hypersurface, then we can naturally extend a family of embeddings by choosing a unit normal $\nu$ and flowing by unit speed, i.e.
$$\varphi_t(x) := \exp_x(t \nu(x))$$
For small enough $t \in (-\epsilon,\epsilon)$, this still gives an embedding $\varphi_t: \Sigma\to M$. Now, one could pull back the metric $g$ by $\varphi$ to obtain a metric on $\Sigma$. This metric gives a (Riemannian) volume form on $\Sigma$, $d\mu_t$. Then,
$$ \frac{\partial}{\partial t} d\mu_t = nH d\mu_t $$
(I think this matches the normalization in @Joseph O'Rourke's answer, but of course the sign changes with the choice direction of the normal)
This follows from the first variation formula.
In other words, the mean curvature measures how the volume of the submanifold (locally) changes under flowing the surface in the direction of the unit normal.
Let $M$ be an oriented hypersurface in an oriented euclidean space $E$.
The normal curvature at a point $x_0\in M$ measures how much it is displaced, in the positive normal direction, with respect to the points of $M$ that are near it, in the same way that the Laplacian of a function measures how much the value of a function at a point differs from the average of the points that are near.
In fact, if $M$ is the graph of a nearly constant function $f:\mathbb R^n\to\mathbb R$ (for example, the vibrating membrane of a drum at some instant) then the Laplacian of $f$ is a good approximation of the normal curvature of $M$.
More generally, if $M$ is possibly not nearly flat, then for any $x_0\in M$ you can find the hyperplane $\Pi$ that is tangent to $M$ at $x_0$, and find a parametrization $f:\Pi\to M$ (only defined on a neighbourhood of $x_0$) such that $p\circ f=id_\Pi$, where $p:E\to\Pi$ is the orthogonal projection. This parametrization shows that $M$ is locally the graph of a function $g=q\circ p$ (where $q$ is the projection on the orthogonal complement of $\Pi$) such that $dg(x_0)=0$. Then, the normal curvature of $M$ at $x_0$ equals $\Delta g(x_0)$.
