It is an actual geophysical problem where we study the liquid flow. We measure at each 2D grid point and each time interval three components of the liquid velocity and we want to compute how the liquid level changes over time.

If I understand the problem correctly (which is beyond this question, unless you specifically know the answer), we can assume that the problem is time-independent, that is we find the liquid level as a function of the spatial coordinates $x,y$ at a specific time $t$.

This problem can be rewritten as an equation for each time step as:

$$ (V_{x}\frac{\partial h}{\partial x}+h\frac{\partial V_{x}}{\partial x}) + (V_{y}\frac{\partial h}{\partial y}+h\frac{\partial V_{y}}{\partial y}) = - V_{z} $$

Here the known quantities are the 3D velocities at each 2D grid point $V_{x},V_{y},V_{z}$.

We can also numerically compute from the known quantities $\frac{\partial V_{x}}{\partial x}, \frac{\partial V_{y}}{\partial y}$.

Our objective is to compute $h=h(x,y)$ at each grid point.

If it was a 1D problem (i.e. only $x$ dependent) then the problem can be solved using the standard numerical techniques derived from the Euler’s method. But how to approach this problem in the 2D case (i.e. $x,y$ dependent)?

To mathematically simplify the problem we can rewrite the equation above like this:

$$ \frac{\partial h}{\partial x}a(x)+hb(x) + \frac{\partial h}{\partial y}c(y)+hd(y) = - e(x,y) $$ where $a,b,c,d,e$ are numerically known and $h=h(x,y)$ is the objective.

  • $\begingroup$ Hmm, it seems that you're converging to a solution, which indicates an iterative method for your matrix equation. $\endgroup$
    – Kyle Kanos
    Mar 19, 2022 at 15:20

1 Answer 1


Given what seems to be the physics of the situation (shallow water incompressible flow with $h$ being the depth) surely we have
$$ V_z(x,y,t)\equiv \frac{\partial h(x,y,t)}{\partial t}. $$ So, given the $V(x,y,t)$'s and $h(x,y,t)$ at at time $t$, you have all the ingredients to find $h(x,y, t+\Delta t)$ by Euler's method.

Added Edit: I see now, after the OP's comments, what is being asked. It is not the time evolution that we want, but to solve have the equation $$ ({\bf v}\cdot \nabla) h(x,y) + h(x,y) (\nabla \cdot {\bf v})= {\rm known} $$ Here ${\bf v}= (v_x, v_y)$ is known, and we have that $h=0$ on some boundary. If we want to find $h$ at the point $(x,y)$ we need to start at $x,y$ and compute (numerically) the curve $x(s), y(s)$ such that $x(0)=x$,$y(0)=y$ and $$ \frac{d x}{ds}= -v_x(x,y)\\ \frac{d y}{ds}= -v_y(x,y). $$ and follow it until at $s=s_0$ it reaches the boundary where $h=0$. These curves--- backtracting the fluid flow lines--- are the characteristics. The we solve the first order equation $$ -\frac{d h}{ds}+ (\nabla\cdot {\bf v})h(s) = v_z(x(s),y(s)) $$ with initial data $h(s_0)=0$. Here the diveregence $(\nabla\cdot {\bf v})$ is known at each $s$ because we can compute $\nabla\cdot {\bf v}$ at each point $x(s) y(s)$ from the given velocity-field data. Then $h(x,y)= h(x(0),y(0))$ solves the problem because the derivative along the curve is $$ \frac{d h}{ds} = -({\bf v}\cdot \nabla) h. $$

Now that we know that the is an answer to the problem we can find $h(x,y)$ much more directly by turning it into a matrix problem. For an $N$-by-$N$ grid Treat the knowns and unknowns as arrays $v_i(n,m)$ ($i=x,y,z$), $h(n,m)$ and the arrays as vectors with $N^2$ entries. Then the problem is linear $$ W h= V_z $$ where $W$ is an $N^2$-by$N^2$ matrix arising from discretizing $$ \nabla\cdot (h {\bf v})= V_z. $$ by nearest neighbour differences. Now invert the matrix to write $h=W^{-1} V_z$. As the $W$ matrix is sparse I would recommend using conjugate-gradient algorithm

  • $\begingroup$ Thank you for your suggestion. However, in my case I do not know h(x,y,t) at at time t=0, I only know 3D velocities at each time t.I can assume that this is a problem with a known boundary condition (h at the boundary is equal to zero), think of it is a viscous fluid for better understanding. The goal is actually to find h(x,y,t) at at time t=0. $\endgroup$
    – SVS
    Mar 22, 2022 at 18:27
  • $\begingroup$ Then you do not have enough enough information to solve the problem. $\endgroup$
    – mike stone
    Mar 22, 2022 at 20:33
  • $\begingroup$ Why do you think so - all parameters in the equation are known, except h(x,y)? I think it is an ordinary linear PDE in relation to x,y. $\endgroup$
    – SVS
    Mar 23, 2022 at 14:52
  • $\begingroup$ Of course! You can use the method of characteristics. You know that $dh/ds=-h(s)( \nabla \cdot {\bf v}_{hor})$ along the characteristic curves $d{\bf r}/ds = {\bf v}_{\rm hor}$. Just integrate the collection of 1-d equations along the curve from each boundary point where $h=0$! Grrr I should have seen that instantly.... $\endgroup$
    – mike stone
    Mar 23, 2022 at 23:51
  • $\begingroup$ Could you please explain a little more your idea? I also feel that I have all required information to solve this problem, and I am also familiar with the analytical method of characteristics, but I cannot find any example on how to use this method in the numerical analysis. I am unfortunately lucking the math background. $\endgroup$
    – SVS
    Mar 24, 2022 at 15:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.