Almost all finite methods that use forward time models adhere to the Courant-Friedrichs-Lewy law which calculates a courant number and compares it to a $C_{max}$, which is what determines stability, so for 2-D:
$$C = \frac {V_xdt}{dx} + \frac{V_ydt}{dy} \ge C_{max}$$
Where $C$ is the courant number, $V_i$ is the velocity in subscripted direction, $dx$ or $dy$ is the length interval in specified direction and $dt$ is the specified time step. If you are getting numerical instability, increasing the mesh spacing will lower the courant number, but then you might be missing the dynamics you want to capture. If you decrease the time step, it will take longer for the dynamics to develop computationally.
One way to solve this computational problem is to code in a dynamic time step: Keep calculating the $C$ and $C_{max}$. If $C > C_{max}$, then have the model lower the step by half. Check again, then reset for the next step so you can have larger timesteps for "less dynamic" times. I have implemented this many times in other models and it has worked well.
Another possibility might be that you have high aspect ratios in your mesh: if your $x$ direction is 10 km and your $y$ is 1 km, a $10:1$ aspect ratio is very high. In your case, since the model package you use are triangles, you want the angles in the triangle grid to be as acute as possible. Obtuse angles would mean high aspect ration, in general.
Like llmari, I am not familiar with hydrodynamic models, so I do not do specifics. My expertise are in both mantle convection and plate flexure models. I hope this helps.