How to simply model fluid flow over heigtmap? My input is a height map, i.e. the ground of the 3D world, represented by a function defined on a $N$-sized grid $h:\{0,N\}^2 \rightarrow \mathbb{R}$. The height of the ground over $(x,y,0)$ is $h(x,y)$.
I would like to model an approximate fluid flow over this height map. As a first attempt, I assume the fluid is perfect, incompressible, and constrained to the ground. A solution to the problem can then be another function $f:\{0,N\}^2 \rightarrow \mathbb{R}$ where $f(x,y)$ is the height of fluid over the ground point $(x,y,h(x,y))$. The fluid top surface is free.
With the perfect and incompressible fluid hypotheses, Navier-Stokes equations are as follow, with $\vec{V}$ the speed vector, $P$ the pressure, and $f$ all other mass forces:
$$\nabla\cdot  \vec{V} = 0$$
$$\frac{\partial\vec{V}}{\partial t} + \vec{V}\cdot\nabla\vec{V}= f - \frac{1}{\rho} \nabla P~~(Eq 1)$$
But I think I can simplify even more these equations due to the last hypothesis, yet I do not know how.

EDIT:
Maybe I could assume
$$\frac {dV_z}{dt} <<< 1$$
and use the projection of the previous equation over $z$ axis to derive, for any $x,y,z$ such as the point $(x,y,z) $ is in the fluid :
$$  P (x,y,z) = P_0 + \rho g (z_{fluid surface} - z) $$
I can then integrate this pressure over the height of the fluid column for a given $x,y $ and find the mean pressure of this column, then use this value in the $x $ and $y $ axes projections of the equation 1 in order to find the speeds.
Would this be a correct reasoning?
 A: I assume you want to solve the shallow water equations, which are the depth-averaged versions of the Euler equations. They can be written as
$\frac{\partial h}{\partial t}+\frac{\partial uh}{\partial x}+\frac{\partial vh}{\partial y}=0$
$\frac{\partial v}{\partial t}+\frac{\partial}{\partial y}\big(v^2+\frac{1}{2}gh^2\big)+\frac{\partial}{\partial x}(uvh)=0$
$\frac{\partial u}{\partial t}+\frac{\partial}{\partial x}\big(u^2+\frac{1}{2}gh^2\big)+\frac{\partial}{\partial y}(uvh)=0$
where $h=h(x,y)$ is the water depth field, $u$ and $v$ are the velocity component fields in $x$ and $y$ directions.
These equations can be solved with many numerical schemes for instance finite difference methods (e.g. using the Lax-Wendroff scheme).
However, your concept of using particles seems to require a Lagrangian approach. From this point of view, the smoothed particle hydrodynamics method (SPH) is a possible choice. Despite of certain difficulties, fortunately, the shallow water equations above are already solved with SPH successfully. I would suggest to read the papers below (if you have access to them) for the technical details of the solution.
[*] Giulia Rossi, Michael Dumbser, Aronne Armanini,
A well-balanced path conservative SPH scheme for nonconservative hyperbolic systems with applications to shallow water and multi-phase flows,
Computers & Fluids,
Volume 154,
2017,
Pages 102-122,
[**] Xilin Xia, Qiuhua Liang, Manuel Pastor, Weilie Zou, Yan-Feng Zhuang,
Balancing the source terms in a SPH model for solving the shallow water equations,
Advances in Water Resources,
Volume 59,
2013,
Pages 25-38
Edit:
Although the implementation of SWE-SPH is easy compared to finite element or finite volume solvers, I think it is better to use an existing solver like SWE-SPHysics, which is an open source tool and a little brother of DualSPHysics. As long as you desire a little more programming and higher configurability, there are further options like Aboria or Nauticle (my own project), being general purpose particle simulation tools.
