Why can't you simulate isotropic fluid flow on a square lattice? There are easy methods for discrete simulations of gas dispersion in two dimensions.  If you take a large square lattice, each cell of which is assumed to contain at most one gas molecule, and you move the molecules from cell to adjacent cell at random, the large-scale results are in many ways a good simulation of gas behavior. In particular, even though molecules are individually moving in only the four cardinal directions, at a large scale the simulation does not show any bias to these four directions. If you start with a large square lump of material, it will rapidly diffuse into a circle.  This property of the simulation is called isotropy.
Long ago when I was involved in research in this area I was told that this worked fine for gas flow simulations, but not for incompressible fluids.  If one wanted to make an isotropic two-dimensional simulation of incompressible fluid flow, I was told, one had to use a hexagonal lattice instead of a square lattice.  I was told that one could deduce from the Navier-Stokes equations that that any simulation of incompressible fluid on a square lattice would necessarily be anisotropic.
Is this correct?  If so, how does the proof go?  If the argument is complex, is it possible to get an intuitive idea of why compressible and incompressible fluids are different in this regard?  What is a reference that would include the full proof?
Addendum: Despite the answer I posted below, I still don't understand any of the details. I would be glad to award the bounty to someone who could explain it.
 A: It seems that this might be the original source: Lattice gas automata for the Navier-Stokes equations, by Frisch et al. (PDF)

The [square lattice] HPP automaton is invariant under $\pi/2$ rotations. Such a lattice symmetry is insufficient to ensure the isotropy of the fourth degree tensor relating momentum flux to quadratic terms in the velocity. …
  Observe that when the underlying microworld is two-dimensional and invariant under the hexagonal rotation group (multiples of $\pi/3$) the tensor $T$ is isotropic…

Unfortunately it will be a long time before I am able to understand this properly; I would still prefer to see an answer that explains what is going on.
A: As you say you are not very familiar with tensors (and I wouldn't call myself an expert either), I will try to explain everything in linear-algebraic language. To quote the paper,

...for the HPP model, the momentum flux tensor has the form
  $$P_{\alpha\beta} = p\delta_{ab} + T_{\alpha\beta\gamma\epsilon}u_\gamma u_\epsilon + O(u^4).$$
  ...The tensor $T$ is, by construction, symmetric in both $(\alpha,\beta)$ and $(\gamma,\epsilon)$.

Here $T$ is "the fourth degree tensor relating momentum flux to quadratic terms in the velocity". We can think of the second-order tensors $P_{\alpha\beta}$ and $u_\gamma u_\epsilon$ as symmetric $2\times2$ matrices ($P_{\alpha\beta}$ must be symmetric to conserve angular momentum, while $u_\gamma u_\epsilon$ is $uu^T$ which is clearly symmetric), and $T$ as a linear transformation mapping the latter to the former.
A tensor is isotropic if it looks the same under any rotation of the coordinate frame. In this case, that means that if $T$ maps some matrix $X$ to a matrix $Y$, and you rotate the coordinate frame so that $X$ becomes $X'$ and $Y$ becomes $Y'$, then $T$ still maps $X'$ to $Y'$. Formally, if a rotation transforms vectors as $v' = Rv$, then matrices transform as $A' = RAR^T$, and we want $T(RXR^T)=RT(X)R^T$.
Things become clearer if we make a judicious choice of basis for the three-dimensional space of symmetric $2\times2$ matrices, namely
$$E_1 = \begin{bmatrix}1 & 0\\0 & 1\end{bmatrix},\quad E_2 = \begin{bmatrix}1 & 0\\0 & -1\end{bmatrix},\quad E_3 = \begin{bmatrix}0 & 1\\1 & 0\end{bmatrix}.$$
Then it turns out that a rotation by an angle $\theta$ transforms this basis as
$$\begin{align}
E_1' &= E_1, \\
E_2' &= E_2\cos2\theta + E_3\sin2\theta, \\
E_3' &= E_3\cos2\theta - E_2\sin2\theta,
\end{align}$$
which one can also see by visualizing them as quadratic functions of the form $x \mapsto x^TEx$ over $\mathbb R^2$. In other words, a rotation transforms $E_2$ and $E_3$ into each other, so for $T$ to be invariant under such a transformation, it has to do "the same thing" to both $E_2$ and $E_3$, in some sense. But if you only look at rotations of $\pi/2$, you don't see this dependence, because then $E_2' = -E_2$ and $E_3' = -E_3$. So a tensor that treated $E_2$ and $E_3$ differently could be invariant under rotations of $\pi/2$, but it would not be invariant under arbitrary rotations. For example, we could choose $T$ such that
$$T(E_1) = T(E_2) = 0,\quad T(E_3) = E_3.$$
You can verify that this tensor is invariant under a rotation of $\pi/2$, but not under a rotation of $\pi/4$ (where $E_2' = E_3$ and $E_3' = -E_2$), so it is not isotropic.

Anyway, since the title of the question still asks "Why can't you simulate isotropic fluid flow on a square lattice?", allow me to reiterate my comment that this is an issue with lattice gas automata in particular, not with fluid simulations on square grids in general. The problem with isotropy comes in when trying to show the convergence of the cellular automaton model to the Navier-Stokes equations. As far as I am aware, you can simulate isotropic fluid flow on a square grid through the usual methods of directly discretizing the Navier-Stokes equations and solving them numerically, without going through an intermediate cellular automata model.
A: I have an existence proof that you can make an asymptotically angularly isotropic three dimensional cellular automata rule that works in a cartesian lattice.  It creates gliders which can move in (asymptotically) any direction at (asymptotically) any velocity. However the rule as written does not deal well with collisions between these gliders. But at least it shows that a fairly simple rule can produce angularly isotropic motion behavior.
https://docs.google.com/document/d/1ocynhGqTrHZ_-GK5EBZjSBAUzj5pPmRlFiYclf1l3-s/edit
