# How to quantify the numerical diffusion term in a second-order upwind advection scheme?

In the first-order upwind scheme, numerical diffusion can be quantified as:

$$\frac{dT}{dt} = -w\left(\frac{dT}{dz}\right) + \left(\frac{wdz}{2}\right)\left(\frac{d^2T}{dz^2}\right)$$

For Lax-Wendroff, it is:

$$\frac{dT}{dt} = -w\left(\frac{dT}{dz}\right) + \left(\frac{w^2dt}{2}\right)\left(\frac{d^2T}{dz^2}\right)$$

It is said that the second order upwind scheme produces less numerical diffusion than the first order counterpart. How can it be quantified for a second-order upwind scheme?

• Can you provide the exact second-order stencil you have in mind? And the precise equation you want to solve? Presumably linear transport? Mar 23 at 16:06
• Hi @Daniel, thanks for answering! I edited the post. I want to calculate the advection of temperature T with ice flow that has a velocity w and grid resolution dz. It is an advection equation (see figures). Mar 23 at 16:18
• Out of curiosity, did you test your scheme? Did it turn out to actually be unstable? Mar 25 at 10:27
• @Daniel Hi, yes I did test it. It is indeed unconditionally unstable. Just like the first order upwind the numerical stability seems to be determined by the CFL criterion. Mar 25 at 16:06

The approach I would try here is that of a modified equation. I am not sure if it goes through for you full problem as your transport speed appears to change with $$k$$ and $$t$$. But I will outline how to analyse diffusion of your scheme when $$w$$ is constant.

First, assuming $$w < 0$$, I rewrite your scheme as

$$T_k^{t+1} = T_k^t - \frac{w \Delta t}{2\Delta z} \left( -T^t_{k+2} + 4 T^t_{k+1} - 3 T^t_{k} \right)$$

Next, to determine the local error, I plug the analytic solution $$T(x,t)$$ into your numerical method. I will play fast and loose with notation but from now on, $$T^t_k$$ stands for $$T(x_k, t)$$. This allows me to do Taylor expansions

$$T^{t+1}_k = T^t_k + \Delta t \partial_t T^t_k + \frac{\Delta t^2}{2} \partial_{tt} T^t_k + \mathcal{O}(\Delta t^3)$$

and

$$T^t_{k+2} = T^t_k + 2 \Delta z \partial_z T^t_k + 2 \Delta z^2 \partial_{zz} T^t_k + \frac{8 \Delta z^3}{6} \partial_{zzz} T^t_k + \mathcal{O}(\Delta z^4)$$

and

$$T^t_{k+1} = T^t_k + \Delta z \partial_z T^t_k + \frac{\Delta z^2}{2} \partial_{zz} T^t_k + \frac{\Delta z^3}{6} \partial_{zzz} T^t_k + \mathcal{O}(\Delta z^4)$$

If I plug these three expansions into the method, a lot of stuff cancels out and I am left with

$$\Delta t \partial_t T^t_k + \frac{\Delta t^2}{2} \partial_{tt} T^t_k + \mathcal{O}(\Delta t^3) = -\frac{w \Delta t}{2 \Delta z} \left( 2 \Delta z \partial_z T^t_k + \frac{2 \Delta z^3}{3} \partial_{zzz} T^t_k + \mathcal{O}(\Delta z^4) \right)$$

Simplifying this gives me

$$\partial_t T^t_k + w \partial_z T^t_k = - \frac{\Delta t}{2} \partial_{tt} T^t_k + \frac{w \Delta z^2}{3} \partial_{zzz} T^t_k + \mathcal{O}(\Delta t^2) + \mathcal{O}(\Delta z^3)$$

Lastly, because $$\partial_{tt} T^t_k = w^2 \partial_{zz} T^t_k$$, we get

$$\partial_t T^t_k + w \partial_z T^t_k = \underbrace{- \frac{\Delta t w^2}{2} \partial_{zz} T^t_k}_{\text{numerical diffusion}} + \frac{w \Delta z^2}{3} \partial_{zzz} T^t_k + \mathcal{O}(\Delta t^2) + \mathcal{O}(\Delta z^3)$$

First, you can see that if your derivatives are all bounded, the error is of order $$\mathcal{O}(\Delta t) + \mathcal{O}(\Delta z^2)$$ which fits a first order explicit Euler time stepper combined with second order finite differences in space.

However, the problem is that the diffusion has the wrong sign: because $$w^2 \Delta t > 0$$, the modified equation for the local error is of the form $$\partial_t T + w \partial_z T = -\alpha \partial_{zz} T$$ with $$\alpha > 0$$ which is unstable.

This indicates that your scheme is unconditionally unstable. Have you tried to implement it to see what it does? I would expect it to blow up after some time, no matter how small you make $$\Delta t$$.