Implementation of 2D advection equation with second order upwind scheme?

Hi all,

I am trying to numerically discretize a 2D advection equation to model the transport of rocks with thickness (h_debris) on top of glacier ice with velocity components (velx_mod and vely_mod). The equation to be discretized is as follows:

dh/dt = -d(u*h)/dx - d(v*h)/dy = -u(dh/dx) - v(dh/dy) - h(du/dx) - h(dv/dy)

% in the model:
% u = velx_mod (x component of ice surface velocity)
% v = vely_mod (y component of ice surface velocity)
% h = h_debris (thickness of the rock layer)

I therefore did the following for a second order upwind scheme in 2D:

if velx_mod(i,j) < 0
    term1_debris_x(i,j) = -velx_mod(i,j).*((-h_debris(i+2,j)+4*h_debris(i+1,j)-3*h_debris(i,j))./(2*deltax_d)) - h_debris(i,j).*((-velx_mod(i+2,j)+4*velx_mod(i+1,j)-3*velx_mod(i,j))./(2*deltax_d));
elseif velx_mod(i,j) >= 0
    term1_debris_x(i,j) = -velx_mod(i,j).*((3*h_debris(i,j)-4*h_debris(i-1,j)+h_debris(i-2,j))./(2*deltax_d)) - h_debris(i,j).*((3*velx_mod(i,j)-4*velx_mod(i-1,j)+velx_mod(i-2,j))./(2*deltax_d));
end

if vely_mod(i,j) < 0
    term1_debris_y(i,j) = -vely_mod(i,j).*((-h_debris(i,j+2)+4*h_debris(i,j+1)-3*h_debris(i,j))./(2*deltax_d)) - h_debris(i,j).*((-vely_mod(i,j+2)+4*vely_mod(i,j+1)-3*vely_mod(i,j))./(2*deltax_d));
elseif vely_mod(i,j) >= 0
    term1_debris_y(i,j) = -vely_mod(i,j).*((3*h_debris(i,j)-4*h_debris(i,j-1)+h_debris(i,j-2))./(2*deltax_d)) - h_debris(i,j).*((3*vely_mod(i,j)-4*vely_mod(i,j-1)+vely_mod(i,j-2))./(2*deltax_d));
end

h_debris(i,j) = h_debrisini(i,j) + deltat_d*(term1_debris_x(i,j) + term1_debris_y(i,j));

I am not an expert in this field and this is the first time I implement this advection equation in 2D. Could someone check my code and assure I did the right thing here? Thanks already.

[Torsten]

u and v are given or are they also solutions of PDEs ?

[Yoni Verhaegen -WE-1718-]

U and V are given, yes!

[Torsten]

If they are given, why do you approximate their derivatives with respect to x and y ?

Can you give an example of their functional form ?

[Yoni Verhaegen -WE-1718-]

The first terms (-u(dh/dx) - v(dh/dy)) consider the advection of thicker or thinner debris from upglacier, the second (-h(du/dx) - h(dv/dy)) represents thickening or thinning of the rock layer due to compressional or extensional ice flow.

[Torsten]

You said u and v are given. How are they given ? As explicit equations depending on x and y ? As discrete data ?

[Yoni Verhaegen -WE-1718-]

The velocity components u and v are produced by a 3D ice dynamic model. I read them in as .nc files to get a velocity field over the entire glacier (see attachment). They are assumed constant over time.

[Torsten]

I have the feeling that keeping velocity constant over time forbids to take into account the terms -h(du/dx) - h(dv/dy) of your equation because thickening or thinning of the rock layer due to compressional or extensional ice flow will itself influence the velocity field. But this was not your question.

So concerning your discretization: the discretization you suggest is a 2nd order upwind method, but maybe not stable since your velocity fields look quite diffuse. I suggest starting with the usual 1st-order upwind scheme. This produces numerical diffusion, but is stable. Further I suggest using (in the case u,v >= 0)

d(u*h)/dx ~ ((u*h)(i)-(u*h)(i-1))/dx instead of u(i)*(h(i)-h(i-1))/dx + h(i)*(u(i)-u(i-1))/dx
d(v*h)/dy ~ ((v*h)(j)-(v*h)(j-1))/dy instead of v(j)*(h(j)-h(j-1))/dy + h(j)*(v(j)-v(j-1))/dy

if you really want to solve

dh/dt + d(u*h)/dx +d(v*h)/dy = 0

instead of

dh/dt + u*dh/dx + v*dh/dy

as I assume to be more correct if u and v are prescribed and kept constant over time (as noted above).

For a stable 2nd order upwind scheme that can cope with discontinuities in u and v see e.g.

https://en.wikipedia.org/wiki/MUSCL_scheme

Furthermore, you will have to think about the boundary conditions for h. For u,v >=0, they have to be prescribed at x=xmin, ymin<=y<=ymax and xmin <=x<=xmax, y=ymin, e.g.

[Yoni Verhaegen -WE-1718-]

Thank you for your input, it really helps me a lot and I appreciate it. So if the velocity field were to change over time, mny original implementation (as stated in the code of my question here), would be correct?

Answer 1

[Torsten]

Formally, your implementation is correct as it is (I did not check the one-sided difference formulae in detail).

The problems I see are:

1. the stability of your scheme (usually, 2nd order schemes which are not stabilized (as your scheme) tend to become instable/oscillatory). That's why you should start with first-order upwind or better: with my suggestion from below.
2. the velocity field. Usually, a velocity field is divergence free to ensure mass-conservation (i.e. the condition du/dx + dv/dy = 0 automatically holds). In this case, your terms h*du/dx + h*dv/dy automatically cancel out and you are left with the equation dh/dt + u*dh/dx + v*dh/dy = 0. You should test this before starting. If you don't get du/dx + dv/dy approximately 0, you must be very careful with the results you get from your simulation.
3. The boundary condition and where you have to set it.
4. The implementation of the scheme. Note that one cell apart from an inflow boundary, you must use a value for h beyond the boundary points (h_debris(i+2,j),h_debris(i-2,j),(-h_debris(i,j+2),h_debris(i,j-2)). So, an extrapolation is necessary.

I strongly recommend using clawpack for your task where a second-order upwind implementation for the 2-dimensional advection equation is already implemented:

https://www.clawpack.org/

[Yoni Verhaegen -WE-1718-]

I got these results with the second order upwind scheme (see attachment), which seems satisfactory to me, no? I am trying to consider all your arguments but for one reason or another, the first order upwind scheme did not produce nice results.

[Torsten]

Looks impressive.

Why does the thickness not decrease in the region from where the growth starts (around x=105, y=72) ? If you integrate your differential equation over the complete region, h should be conserved (i.e. shouldn't become "more" or "less" in the course of time, only change its position).

[Yoni Verhaegen -WE-1718-]

It is because for this example, I assumed a steady deposition rate, i.e. the same amount of rock is repeated on the input pixel every timestep instead of just 1 deposit at the first timestep.

[Torsten]

Then your equation must include a source term s in its formulation:

dh/dt = -d(u*h)/dx - d(v*h)/dy + s(x,y)

And did you check the validity of the global balance :

integral integral h(x,y,t) dx dy - integral integral h(x,y,0) dx dy = t* integral integral s(x,y) dx dy

?

[Yoni Verhaegen -WE-1718-]

Unfortunately, I am not an expert in this field, so I don't really know how I neede to check and handle these issues. I put a lot of effort in it already to get it where it is now, and for that I am really happy that you have been helping me. Would it be an idea to share the model code somehow so that we can figure it out together?