Why does the volume element in spherical polar coordinates contain a sine of the zenith angle?

Since you (the OP) haven't accepted an answer, I'm posting this, but consider this as a supplement to amd's answer, since his/her contribution made me understood this problem, about which I was recurrently thinking for two days. Also the comment of Jahan Claes helped, which was the same that one of my teachers told me. I would also like to note, that the reasoning given here is an intuitive, physicist's reasoning, perhaps not really fit for a mathematics site. However, I think that the question also refers to acquiring the volume element with this kind of reasoning, in an intuitive way.

Consider a point in 3D space described by coordinates $(r,\theta,\phi)$ such as in Fig. 1. Figure 1.

Now consider an infinitesimal increase $dr$ in $r$. Whatever the value of $theta$ and $phi$ are, the length of this increase in the geometrical space will be the same. See Fig. 2. Figure 2.

Now, let's forget about the increase in $r$ for a moment, and consider an infinitesimal increase in $\theta$, and let's just foucus on that. Whatever the value of $\phi$ is, it won't change the arc length, but it strongly depends on the value of $r$, namely as $rd\theta$. You can see this by intersecting the 3D space with a plane perpendicular to the $xy$ plane, and going through the line which connets the orgin $(0,0,0)$ and the point $(r,\theta,\phi)$. For illustration see Fig. 3.

Figure 3.

Let's forget about the increse in $\theta$ too, and let's just focus on what happens, if there is and infinitesimal increase $d\phi$ in $\phi$. To see what happens, I suggest looking on the setup from "above", as shown in Fig. 4. (the $z$ axis is pointing outside of the screen). You can see now, that if $\theta$ is smaller ($\theta_1$), the arclength due to the increase in $\phi$ will also be smaller. However, if $\theta$ is larger, ($\theta_2$), than from above you can see that the arclength will also be larger. If you determine the exact relation, you will see that the arc length is proportional to the projection of the position vector on the $xy$ plane, and $d\phi$. That is the arclength is $|Proj_{\theta}(\vec{r})|d\phi=r\sin\theta d\phi$. Figure 4.

Now, to put it all together, refer to Fig. 5. The volume element comes from the infinitesimal increase in $r$, $\phi$, and $\theta$. This is of course a trapesoidlike thing with pieces of spherical surfaces on its top and bottom. However, since it is infinitely small, it may be approximated by a cuboid (the precise justification of this approximation I can't do for the moment), and its volume may be written as the length of its sides multiplied by eachother, that is its volume is \begin{equation}dV=drrd\theta r\sin\theta d\phi=r^2\sin\theta drd\theta d\phi \end{equation} Figure 5.

Two sidenotes:

1. I really like amd's idea, the same one that one of my teachers told me. That is the term $r^2\sin\theta$ which we have acquired intuitively, is the so called "volume distortion factor", which shows you that if you were to parametrize the geometrical space with three numbers denoted by $r$, $\theta$ and $\phi$, than the parameter space would have three axes denoted by the same letters. If in that paremeter space you would take an infinitesimal cuboid with sides $dr$, $d\theta$ and $d\phi$, with one of its corners residing at $(r,\theta,\phi)$, than the volume of this infinitesimal increment after transforming it back to the real, geometrical space would appear as $r^2\sin\theta drd\theta d\phi$. In genral this volume distortion factor is given by the Jacobian. For the case of transforming spherical coordinates into Descartes coordinates, the value of the Jacobian is indeed $r^2\sin\theta$.
2. One may argue, that in my approximation of the strange trapesoid-like shape as a cuboid, I could have taken the length of the upper arches in computing its volume, that is its volume could have been $dV=dr*(r+dr)d\theta*(r+dr)sin\theta d\phi$. Or one could also argue why I have used $\theta$ in the calculations, not $\theta+d\theta$. That is the volume could also have been $dV=dr*(r+dr)d\theta*(r+dr)sin(\theta+d\theta) d\phi$. However the limiting process done while integrating makes those terms disappear, which contain the squares of infinitesimally small elements, if they appear in a sum that contains other infinitesimally small elements on a lower power. That is e.g. $const*dr+const*dr^2\approx const*dr$, $const*r+const*dr\approx const*r$, and so on. In the second case of computing the volume element, where I have used the infinitesimally increased $\theta$, one would also have to use a trigonometric identity on the sine, and use the approximations $\sin d\theta\approx d\theta$ and $\cos d\theta\approx 1$ to see this.

It's because all the longitudes of a sphere converge toward the north (and south) pole.

The length $L$ in the picture falls as $\sin(\theta)$ when you move up toward the north pole, i.e. when you decrease $\theta$ from $\pi /2$ to $0$ , you can prove this geometrically by finding the relationship between $L=r\mathrm{d}\phi \times\sin(\frac{\pi}{2})$ and the new length of $L\,$ for a given $\theta <\pi /2$:

The angle $\mathrm{d}\phi$ remains constant as we decrease $\theta$.

The volume element is simply $\,r\sin(\theta)\mathrm{d}\phi\times r\mathrm{d}\theta \times\mathrm{d}r$ .

Here’s a way of looking at the coordinate transformation that might help you develop an intuition for this. Instead of thinking of spherical coordinates as a different way of labeling the points in a Euclidean space, think of the spherical-to-Cartesian transformation as defining a mapping between two Euclidean spaces: one in which the coordinate axes are labeled $r$, $\phi$, $\theta$ and in the other $x$, $y$, $z$. From this point of view, when you change an integral over some region $W$ in the $xyz$-space to spherical coordinates, you transform it into an integral over a region $R$ in $r\phi\theta$-space, where the image of $R$ under the mapping is $W$. Now, the volume element in $r\phi\theta$-space is just $dr\;d\phi\;d\theta$, but for the transformed integral you have to use the image of this volume element, so there’s going to be a scale factor involved.

To see what this scale factor might be, examine what happens to the image of a small parallelepiped as you move it around in $r\phi\theta$-space. As you move away from the origin parallel to the $r$-axis, the image of that small volume moves away from the $xyz$-space origin and grows. That’s where a factor of $r$ comes in. If you vary $\phi$, neither the distance of the image from the $xyz$-space origin nor the size of the image change, thus $\phi$ doesn’t appear in the scale factor. If you vary $\theta$, the image of the small volume stays the same distance from the $xyz$-space origin, but moves closer to and farther away from the $z$-axis, shrinking and growing as it does. The size of the image is proportional to its distance from the $z$-axis, but that distance isn’t a linear function of $\theta$. In the diagram in Aniket’s answer, this distance is labeled $\rho$, which is obviously equal to $r\sin\theta$.