Hello,
Looking again at your case, I noticed that you have several overconstrained triangles in the mesh, for example in the North-West corner. It means that the 3 points of a triangle are on a border. Normally the mesh-generators should avoid this (Matisse for example had an option for suppressing such triangles). You can also rework on the mesh and split them into two triangles. The problem is that some numerical schemes are sensitive to this (and others absolutely not...). Just to make you understand, imagine that the boundary conditions impose the velocities of the 3 points of the triangle, and that you want to solve div(U)=0, this may be impossible for this triangle, as all the velocities are already imposed, which leaves you no degree of freedom.
This may explain that the diffusion step has badly conditioned linear systems. We have actually many other cases with tidal flats that raise no problem with linear systems.
With best regards,
Jean-Michel Hervouet