I motivate these grad, curl, and div for myself as being the things which would make the respective version of Stokes' theorem true infinitesimally. You can read about this interpretation in the context of differential forms here:
https://math.stackexchange.com/a/614473/34287
Given that interpretation, we can see
- The gradient at a point $p$ is the vector which represents the map $df|_p$. In other words $df|_p(\vec{v}) = \nabla f|_p \cdot \vec{v}$. You can see that the dot product is maximized when $\vec{v} = \nabla f|_p$, which justifies the "steepest slope" interpretation.
- Given a vector field $F$, convert it into a covector field $\omega$ via $\omega(\vec{v}) = F \cdot \vec{v}$. Then the integral of the differential form $\omega$ along a curve is the same as the "line integral" of the vector field $F$. Following my answer above, you can define $d\omega$ in terms of the line integral around little parallelograms. $d\omega$ is a 2-form, but it can be converted into a vector field by a similar representation trick (this one only works in $\mathbb{R}^3$): define $\textrm{curl}(F)$ as the vector for which $\textrm{curl}(F) \cdot (\vec{v} \times \vec{w}) = d\omega(\vec{v},\vec{w})$. The RHS is maximized when the line integral around a small parallelogram defined by $\vec{v}$ and $\vec{w}$ is maximized and then the LHS would require that $\textrm{curl}(F)$ be perpendicular to that parallelogram, and pointing in the direction given by the right hand rule so as to make the integral positive.
- Given a vector field $F$, convert it into a $2$-form $\omega$ via $\omega(\vec{v},\vec{w}) = \textrm{Det}(\vec{v},\vec{w},\vec{F}) = \vec{F} \cdot ( \vec{v} \times \vec{w})$. The integral of $\omega$ over a surface is the same as the surface integral of $F$. We can define $d\omega(\vec{v},\vec{w},\vec{b})$ in terms of the integral of $\omega$ around the boundary of the small oriented parallelepiped defined by the three vectors. This is the same as the surface integral of $F$ over that same surface. This is a top level form in $\mathbb{R}^3$, so it is a multiple of the volume form, i.e. $\textrm{div}(F)\textrm{Det}(\vec{v},\vec{w},\vec{b}) = d\omega(\vec{v},\vec{w},\vec{b})$. So the sign of the divergence can be "seen" according to whether most of the vectors $F$ are pointing "into" or "out of" a small parallelepiped based at $p$.
You can replace the parallelogram with a circle and the parallelepiped with a sphere to obtain a more "symmetric" representation, but you loose the easy connection with the differential forms. Also, it is much easier to decompose surfaces into small parallelograms (and solids into small parallelepipeds), which leads to the global Stokes' theorems being just "telescoping" sums, with cancellation of all interior terms.
I think the real difficulty with motivating these things comes from the "unnatural" conversion of differential forms into scalar and vector fields, in ways which only work in $\mathbb{R}^3$. They are somewhat arbitrary. They are also unnecessary, as any work you want to do can be done with the differential forms themselves.
These will let you draw pictures though! In each case, you have to imagine many different "test" domains of integration, and try to figure out which domain will maximize the integral. The grad, curl, and div connect to those maximizing domains.
EDIT: I forgot to mention the laplacian, which doesn't really fit into this framework. The interpretation you link to seems pretty much optimal though.