What's the intuitive reason that phase space flow is incompressible in Classical Mechanics but compressible in Quantum Mechanics? One of the most important results of Classical Mechanics is Liouville's theorem, which tells us that the flow in phase space is like an incompressible fluid. 
However, in the phase space formulation of Quantum Mechanics, one of the main results due to Moyal is that quantum flows are compressible.
So what's the intuitive reason for this difference?
Formulated a bit differently: What exactly is the assumption used in the derivation of Liouville's theorem that is no longer valid in Quantum Mechanics?
 A: *

*The so-called (generic) failure of the quantum Liouville theorem, i.e. the (generic) violation of the continuity equation
$$ \rho~{\rm div}_{\rho} X^Q_{-H} + \frac{\partial \rho}{\partial t}~\neq~0 \tag{1}$$
$$ \begin{align}\Leftrightarrow\qquad \frac{\partial \rho}{\partial t}
~\stackrel{(1)}{\neq}~&\rho~{\rm div}_{\rho} X^Q_H
~\stackrel{(6)+(7)}{=}~\rho~{\rm div}_{\rho} X_H \cr
~\stackrel{\text{Leibniz}}{=}& X_H[z^I] ~\frac{\partial\rho}{\partial z^I}
~+~ \rho~\underbrace{ {\rm div}_1 X_H}_{=0}\end{align}\tag{2} $$
for the quantum flow on a $2n$-dimensional phase space, can intuitively be understood as the appearance of higher-order differential operators $X$ in the $\star$-product, which do not (necessarily) obey Leibniz rule
$$ X[fg] ~\neq~X[f]g+f X[g] .\tag{3}$$ 

*In eq. (1) we have defined the quantum ($Q$) version 
$$ X^Q_H~:=~\frac{1}{i\hbar}[H\stackrel{\star}{,}\cdot] ~=~
\frac{2}{i\hbar}H\sinh\left(\frac{i\hbar}{2}\stackrel{\leftarrow}{\partial} \wedge\stackrel{\rightarrow}{\partial}\right) ~=~X_H + {\cal O}(\partial^3) \tag{4}  $$ 
of a Hamiltonian vector field
$$ X_H ~:=~ \{H, \cdot \}~=~H\stackrel{\leftarrow}{\partial} \wedge\stackrel{\rightarrow}{\partial} .\tag{5}$$ 
Note that the coordinate components are the same
$$ X^Q_H[z^I]~\stackrel{(4)}{=}~X_H[z^I],\tag{6}$$
which is part of the problem. Also in eq. (1) we have for convenience defined a divergence 
$${\rm div}_{\rho} X~:=~ \rho^{-1}\frac{\partial(\rho X[z^I])}{\partial z^I}\tag{7}$$
of a possibly higher-order differential operator $X$. Eq. (7) is not a geometric object, which foretells the doom of eq. (1).  

*The quantum Liouville theorem (1) is replaced by the quantum Liouville equation
$$ 0~=~\frac{d\rho}{dt}~=~X^Q_{-H}[\rho]+\frac{\partial \rho}{\partial t}\tag{8}$$
$$ \qquad\Leftrightarrow\qquad\frac{\partial \rho}{\partial t}
~\stackrel{(8)}{=}~X^Q_{H}[\rho]
~\stackrel{(3)}{\neq}~ X^Q_H[z^I] ~\frac{\partial\rho}{\partial z^I} ~\stackrel{(6)}{=}~ X_H[z^I] ~\frac{\partial\rho}{\partial z^I}  .\tag{9} $$
The ineq. (9) is precisely the ineq. (2). $\Box$

*See also this related Phys.SE post.
A: Apologies for my inability to share intuition, a frequently subjective issue... I have learned a lot by reading the Steuernagel group numerical flows and topological features of such flows, in practice.  For a recent discussion/proof of the zeros, singularities,and  negative probability density features, hence your source-sink query in anharmonic quantum systems see Kakofengitis, Oliva & Steuernagel, 2017. Basically all bets are off when you (a point in phase space) and  neighbors enter a phase-space cell of order $\hbar$, by dint of the uncertainty principle, and that includes a definition of what a trajectory is. 
If you watch the nifty movies of Cabrera and Bondar in the WP article you are linking to for the Morse and quartic potentials, you actually see this in real time, as a lump (you) spread all over the phase space in a highly organized way... I defy you to discern trajectories there! There is powerful topology at work, but I'd defer to Steuernagel for that.
As a practical reassurance, I'll work out a trivial exercise from our book, on compressibility of Euler flows.  For a Hamiltonian $H=p^2/(2m)+V(x)$, the  Moyal evolution equation amounts to an Eulerian probability transport continuity equation, 
$$\frac{\partial f(x,p)}{\partial t} +\partial_x J_x + \partial_p J_p=0~,$$
where, for $\mathrm{sinc}(z)\equiv \sin z/~ z$ , the phase-space flux is 
$$J_x=pf/m~ ,\\ 
J_p= -f \mathrm{sinc} \left( {\hbar \over 2} \overleftarrow {\partial _p} \overrightarrow {\partial _x} \right)~~ 
\partial_x V(x).
$$
Classical mechanics is crucially different, in that the phase-space current is always $ {\bf J}=(p/m,-\partial_xV(x))f$, and the velocity ${\bf v}=(p/m,-\partial_xV(x))$, manifestly divergenceless in phase space. 
Now note for the oscillator, $V_1=  x^2/2$, $J_p=-f   x  $, so the phase-space velocity ${\bf v}=(mp,- x)$ and $\nabla \cdot {\bf v}=0$, incompressibility. This is a reminder that the quantum oscillator is basically classical, and its wave packets do not spread, as iconically pointed out by Schroedinger... coherent states.  But this is a crying exception.
For a more generic potential, like the quartic, $V_2=  x^4/4$, 
$$
v_p=J_p/f=   -x^3 +\hbar^2 x ~\partial_p^2f /f , \\
\nabla \cdot {\bf v}=  \hbar^2 x ~\partial_p (\partial_p^2 f(x,p) /f(x,p))\neq 0,
$$ 
so the flow is modified by $O(\hbar^2)$ to compressible.
So the strictly quantum difference between the quantum Moyal bracket and the classical Poisson bracket is the crucial element in increasing or decreasing the amount of (quasi)probability in a comoving phase-space region $\Omega$, since
$$
{d  \over dt}\! \int_{\Omega}\! \! dx dp ~f= 
\int_{\Omega}\!\!   dx dp \left ({\partial f \over \partial t}+
  \partial_x (\dot{x} f) + \partial_p (\dot{p} f ) \right)
 =   
\int_{\Omega}\! \!\! dx dp~ (\{\!\!\{ H,f\}\!\!\}-\{H,f\})\neq 0 ~. 
$$


*

*Note added: It is even odder. Quantum flows display physically significant viscosity!

A: Here's my really naive attempt to answer my own question. Please correct me where I'm wrong.
Each point in phase space corresponds to one specific state of the system $(\vec{q}_1,\vec{q}_2,\ldots,\vec{p}_1,\vec{p}_2,\ldots)$. As time passes on this point moves and traces out an orbit in phase space. This orbit can be calculated using Hamilton's equations.
Neighboring points describe similar states. So when we are uncertain about the exact state of our system (which we always are thanks to our limited measurement accuracy) we have to take this into account by using a phase space distribution function. This function $\rho (p,q)$ determines the probability $\rho (\vec{q}_1,\vec{q}_2,\ldots,\vec{p}_1,\vec{p}_2,\ldots)\,d^{n}q\,d^{n}p$ that the system will be found in the infinitesimal phase space volume $d^{n}q\,d^{n}p$. Liouville's equation determines the "orbit" of our initial phase space distribution function. The path that is traced out this way defines a flow in phase space. The two essential ingredients in the derivation of Liouville's equation are 


*

*Hamilton's equations

*The continuity equation for $\rho (p,q)$.


The second ingredient here leads us to the famous conclusion that the phase space flow is incompressible. What this means is that we can put down a pencil on each possible initial configuration (possible in a statistical sense since we are not 100% sure about the initial configuration) and then trace out the phase space flow by moving these pencils through our phase space.
Now in Quantum Mechanics, this is no longer true. Our phase space flow is compressible. In other words, the continuity equation is no longer correct since there are sources and sinks.  What this means is that when we try to trace our obits using pencils we will fail. One trajectory can split into two and other trajectories possibly vanish. (There are sources for new trajectories and sinks where trajectories end.) 
This is a result of the fundamental uncertainty in Quantum Mechanics. While there can also be uncertainty in Classical Mechanics (which is why we use the probability function and Liouville's equation in the first place) it is of a different kind. In Quantum Mechanics there isn't one unique orbit for each possible initial configuration. This is what we mean when we say the phase space flow in Quantum Mechanics is compressible. 
