4.3 Spatial Advection Schemes

Learning Objectives

After completing this lecture, students should be able to:

understand how different numerical schemes approximate spatial derivatives
implement upwind, central-difference, and spectral advection schemes
explain the role of cyclic boundary conditions in numerical advection
compare numerical diffusion and numerical dispersion
evaluate the strengths and weaknesses of different advection schemes
select an appropriate advection scheme for a given application

Concept

Why Different Advection Schemes?

The one-dimensional advection equation is

\[ \frac{\partial T}{\partial t} = F_{\rm adv}, \]

where

\[ F_{\rm adv} = -u \frac{\partial T}{\partial x}. \]

The key numerical challenge is computing the spatial derivative,

\[ \frac{\partial T}{\partial x}. \]

Different numerical schemes approximate this derivative in different ways.

The choice of scheme affects:

accuracy
stability
computational cost
numerical diffusion
numerical dispersion

Role of Cyclic Boundary Conditions

In this section, we use a cyclic boundary condition so that information leaving one side of the domain re-enters from the opposite side.

This choice is useful because all advection schemes require neighboring information near the boundaries.

The upwind scheme requires information from the upstream side of the flow.

The central-difference scheme requires information from both sides of each grid point.

The spectral scheme assumes a periodic domain because Fourier basis functions are cyclic.

Therefore, cyclic boundaries provide a consistent framework for comparing different advection schemes.

Upwind Scheme

The upwind scheme uses information from the upstream direction of the flow.

For positive velocity (\(u>0\)),

\[ \frac{\partial T}{\partial x} \approx \frac{T_i-T_{i-1}} {\Delta x}. \]

Advantages:

simple
stable
robust

Disadvantages:

excessive numerical diffusion
smoothing of sharp gradients

As a result, the bell-shaped tracer gradually broadens and weakens during advection.

Central-Difference Scheme

The central-difference scheme uses information from both sides of a grid point:

\[ \frac{\partial T}{\partial x} \approx \frac{T_{i+1}-T_{i-1}} {2\Delta x}. \]

Advantages:

higher accuracy
reduced numerical diffusion

Disadvantages:

can generate oscillations
may become unstable

Because the scheme does not contain built-in numerical damping, small numerical errors may grow with time.

Spectral Scheme

Instead of computing derivatives using neighboring grid points, spectral methods represent the solution as a sum of waves.

For a periodic domain,

\[ T(x) = \sum_k \hat{T}_k e^{ikx}. \]

The spatial derivative becomes

\[ \frac{\partial T}{\partial x} = \sum_k ik\hat{T}_k e^{ikx}. \]

The derivative is computed exactly in spectral space and then transformed back to physical space.

Advantages:

very high accuracy
minimal numerical diffusion
ideal for smooth solutions

Disadvantages:

requires periodic boundaries
more complicated implementation
less effective for discontinuities

Many global atmospheric and oceanic models use spectral methods because of their excellent accuracy.

Numerical Diffusion and Numerical Dispersion

Numerical schemes introduce errors that can modify the solution.

Numerical diffusion:

smooths gradients
broadens peaks
reduces amplitude

Numerical dispersion:

distorts wave propagation
produces oscillatory noise
shifts phase speeds

Upwind schemes tend to produce numerical diffusion.

Central-difference schemes tend to produce numerical dispersion.

Spectral methods generally minimize both errors for smooth solutions.

Model / Equation

Upwind Scheme with Cyclic Boundary

For positive velocity (\(u>0\)),

\[ F_{\rm adv} = -u \frac{T_i-T_{i-1}} {\Delta x}. \]

At the left boundary,

\[ T_{-1} = T_{N-1}. \]

In Python:

dTdt = np.zeros_like(T)

dTdt[1:] = -U * (T[1:] - T[:-1]) / DX
dTdt[0]  = -U * (T[0] - T[-1]) / DX

Central-Difference Scheme with Cyclic Boundary

The central-difference approximation is

\[ F_{\rm adv} = -u \frac{T_{i+1}-T_{i-1}} {2\Delta x}. \]

This scheme requires information from both sides of each grid point. Therefore, two boundary values must be specified: one for the left side and one for the right side.

For a cyclic boundary, values outside the domain are supplied from the opposite side of the domain:

\[ T_{-1} = T_{N-1}, \]

\[ T_N = T_0. \]

In Python, the interior points can be computed using array slicing:

dTdt = np.zeros_like(T)

dTdt[1:NX-1] = -U * (T[2:NX] - T[0:NX-2]) / (2.0 * DX)

Then the two boundary points are computed using the cyclic boundary condition:

dTdt[0] = -U * (T[1] - T[-1]) / (2.0 * DX)
dTdt[-1] = -U * (T[0] - T[-2]) / (2.0 * DX)

Here, T[-1] provides the left-side neighbor for T[0], and T[0] provides the right-side neighbor for T[-1].

Combining the interior and boundary calculations gives:

def advection_tendency(T, x, U):

    dTdt = np.zeros_like(T)
    NX = len(T)

    dTdt[1:NX-1] = -U * (T[2:NX] - T[0:NX-2]) / (2.0 * DX)

    dTdt[0] = -U * (T[1] - T[-1]) / (2.0 * DX)
    dTdt[-1] = -U * (T[0] - T[-2]) / (2.0 * DX)

    return dTdt

This implementation avoids np.roll() and explicitly shows how the cyclic boundary condition supplies the missing neighboring values.

Spectral Scheme with Cyclic Boundary

The tracer field is represented as a sum of Fourier waves:

\[ T(x) = \sum_k \hat{T}_k e^{ikx}. \]

The derivative in spectral space is

\[ \widehat{ \frac{\partial T}{\partial x} } = ik\hat{T}_k. \]

The advection tendency becomes

\[ F_{\rm adv} = -u \mathcal{F}^{-1} \left( ik\hat{T}_k \right). \]

In Python:

Tk = np.fft.fft(T)

k = 2.0 * np.pi * np.fft.fftfreq(len(T), d=DX)

dTdx = np.real(
    np.fft.ifft(
        1j * k * Tk
    )
)

dTdt = -U * dTdx

Because Fourier basis functions are periodic, the cyclic boundary condition is built directly into the spectral derivative.

Sample Code

The template file for this section is

Template: 4_3-advection3_template.py

Run the script:

python 4_3-advection3_template.py

The script produces

4_3-advection3_template.mp4
4_3-advection3_template.png

This template uses the one-dimensional advection equation with a cyclic boundary condition. The default numerical scheme is the upwind scheme.

The key part of the template is the advection_tendency() function:

def advection_tendency(T, x, U):

    dTdt = np.zeros_like(T)

    dTdt[1:] = -U * np.diff(T) / DX

    # Cyclic Boundary
    dTdt[0] = -U * (T[0] - T[-1]) / DX

    return dTdt

In the exercises, you will replace this function with different spatial advection schemes and compare the results.

Exercise

Step 1. Run the Upwind Scheme

First, create a copy of the template script:

cp 4_3-advection3_template.py 4_3-advection3_upwind.py

Run the model:

python 4_3-advection3_upwind.py

Examine the output files:

4_3-advection3_upwind.mp4
4_3-advection3_upwind.png

The default configuration uses

NT = 200
DT = 1.0
DX = 1.0
U = 0.5

The cyclic boundary condition allows the bell-shaped anomaly to leave the right side of the domain and re-enter from the left side.

Questions:

Does the anomaly propagate toward positive \(x\)?
What happens when the anomaly reaches the right boundary?
Does the anomaly re-enter from the left boundary?
Does the amplitude decrease with time?
Does the shape become smoother or broader?

This case will be used as the reference simulation.

Step 2. Implement the Central-Difference Scheme

Next, replace the upwind scheme with the central-difference scheme.

Create a new script:

cp 4_3-advection3_upwind.py 4_3-advection3_central.py

Replace the advection_tendency() function with:

def advection_tendency(T, x, U):

    dTdt = np.zeros_like(T)
    NX = len(T)

    # Interior points
    dTdt[1:NX-1] = -U * (T[2:NX] - T[0:NX-2]) / (2.0 * DX)

    # Cyclic boundary points
    dTdt[0] = -U * (T[1] - T[-1]) / (2.0 * DX)
    dTdt[-1] = -U * (T[0] - T[-2]) / (2.0 * DX)

    return dTdt

Run the model with the default time step:

DT = 1.0
NT = 200

python 4_3-advection3_central.py

Questions:

Does the central-difference solution remain stable?
Do oscillations appear?
Does the amplitude grow with time?
How does this behavior differ from the upwind scheme?

Step 3. Reduce the Time Step

The central-difference scheme produces less numerical diffusion than the upwind scheme, but it can become numerically unstable when combined with the Forward Euler time integration scheme.

One possible approach is to reduce the time step.

First, use:

DT = 1.0
NT = 200

Run the model and examine the solution.

Next, reduce the time step:

DT = 0.5
NT = 200

Run the model again.

Questions:

Does the instability become weaker?
Does the simulation cover the same total physical time interval?
How does the solution compare with the DT = 1.0 case?

Because the total integration time is now shorter, increase the number of time steps:

DT = 0.5
NT = 400

Run the model again.

Questions:

Does the instability still appear?
Is the instability delayed compared with the DT = 1.0 case?
Does reducing the time step completely eliminate the instability?

Finally, try a smaller time step:

DT = 0.2
NT = 1000

Run the model again.

Questions:

Does further reducing DT improve the solution?
Is the solution completely stable?
What does this suggest about the Forward Euler time integration scheme?

This experiment demonstrates that reducing the time step may improve the solution and delay instability, but it does not fundamentally solve the stability problem associated with the Forward Euler method.

Step 4. Improve Time Integration Using RK4

Although reducing DT improves the numerical solution, it requires substantially more time steps and increases computational cost.

Instead of simply decreasing the time step, we can improve the time integration method.

Create a new script:

cp 4_3-advection3_central.py 4_3-advection3_central_rk4.py

Replace the Forward Euler update:

dTdt = advection_tendency(T, x, U)

T = T + DT * dTdt

with the fourth-order Runge-Kutta (RK4) method:

k1 = advection_tendency(T, x, U)
k2 = advection_tendency( T + 0.5 * DT * k1, x, U)
k3 = advection_tendency( T + 0.5 * DT * k2, x, U)
k4 = advection_tendency( T + DT * k3, x, U)

T = T + DT / 6.0 * ( k1 + 2.0 * k2 + 2.0 * k3 + k4)

Restore the original time step:

DT = 1.0
NT = 200

Run the model:

python 4_3-advection3_central_rk4.py

Questions:

Does RK4 improve the stability of the central-difference scheme?
Does RK4 reduce the need for a very small DT?
How does the RK4 solution compare with the DT = 0.2 Forward Euler solution?
Does RK4 preserve the shape of the tracer better than the upwind scheme?

This experiment demonstrates that numerical stability depends on both the spatial discretization and the time integration method.

Step 5. Implement the Spectral Scheme

Next, replace the central-difference derivative with a spectral derivative while retaining the RK4 time integration.

Create a new script:

cp 4_3-advection3_central_rk4.py 4_3-advection3_spectral.py

Replace the advection_tendency() function with:

def advection_tendency(T, x, U):

    NX = len(x)
    L = x[-1] - x[0]

    T_hat = np.fft.fft(T)

    if NX % 2 == 0:
        k = np.asarray( list(range(0, NX // 2)) + [0] + list(range(-NX // 2 + 1, 0)), dtype=np.float64)
    else:
        k = np.asarray( list(range(0, (NX - 1) // 2)) + [0] + list(range(-(NX - 1) // 2, 0)), dtype=np.float64)

    k = k * 2.0 * np.pi / L
    dTdx = np.fft.ifft( 1.0j * k * T_hat).real
    dTdt = -U * dTdx

    return dTdt

Run the model:

python 4_3-advection3_spectral.py

Examine the output files:

4_3-advection3_spectral.mp4
4_3-advection3_spectral.png

Questions:

How does the spectral solution compare with the upwind solution?
How does it compare with the central-difference RK4 solution?
Which scheme best preserves the shape of the bell-shaped anomaly?
Does the spectral method introduce numerical diffusion?
Why does the spectral method require a cyclic boundary condition?

Step 6. Compare All Schemes

Compare the output animations and final figures from:

4_3-advection3_upwind.mp4
4_3-advection3_central.mp4
4_3-advection3_central_rk4.mp4
4_3-advection3_spectral.mp4

and

4_3-advection3_upwind.png
4_3-advection3_central.png
4_3-advection3_central_rk4.png
4_3-advection3_spectral.png

Focus on the following comparisons:

Scheme	Spatial Scheme	Time Scheme
Upwind	Upwind	Forward Euler
Central	Central Difference	Forward Euler
Central RK4	Central Difference	RK4
Spectral	Spectral	RK4

Questions:

Which scheme produces the strongest numerical diffusion?
Which scheme best preserves the amplitude?
Which scheme best preserves the shape?
Which scheme produces oscillatory or unstable behavior?
How much improvement does RK4 provide relative to Forward Euler?
Which scheme appears most accurate for this smooth periodic tracer?

This comparison highlights the importance of both spatial discretization and time integration in numerical modeling.

Optional: Summary Code

Students who would like a cleaner implementation may use:

4_3-advection3_summary.py

This script combines:

upwind, central-difference, and spectral advection schemes,
Forward Euler and RK4 time integration,
cyclic boundary conditions,

within a single framework.

The user can switch among schemes by modifying a small number of lines in the code. This summary version will also serve as a useful starting point for later chapters when additional physical processes are introduced.

Key Points

The advection tendency

\[ F_{\rm adv} = -u \frac{\partial T}{\partial x} \]

can be approximated using several numerical schemes.

The upwind scheme uses information from the upstream direction of the flow and is generally stable but introduces numerical diffusion.
Numerical diffusion smooths gradients, broadens peaks, and reduces amplitudes.
The central-difference scheme is more accurate than the upwind scheme because it uses information from both sides of a grid point.
The central-difference scheme contains much less numerical diffusion but can become unstable when combined with Forward Euler time integration.
Reducing the time step can delay numerical instability, but it may not completely eliminate the problem.
Numerical accuracy depends on both the spatial discretization and the time integration scheme.
The RK4 method improves numerical stability and accuracy by evaluating the tendency four times during each time step.
Spectral methods compute spatial derivatives using Fourier transforms rather than finite differences.
Spectral methods are highly accurate for smooth periodic fields and require cyclic boundary conditions.
For smooth advection problems, the combination of spectral derivatives and RK4 time integration often provides the most accurate solution.
Numerical models must balance accuracy, stability, and computational cost.

Homework

Complete the following tasks.

1. Central-Difference Scheme with RK4

Use the central-difference spatial derivative together with the RK4 time integration scheme.

Generate and submit:

4_3-advection3_central_rk4.mp4

and

4_3-advection3_central_rk4.png

In a short paragraph, describe:

how the solution differs from the Forward Euler central-difference solution,
whether RK4 improves numerical stability,
whether numerical diffusion is reduced relative to the upwind scheme.

2. Spectral Scheme with RK4

Use the spectral derivative together with the RK4 time integration scheme.

Generate and submit:

4_3-advection3_spectral.mp4

and

4_3-advection3_spectral.png

In a short paragraph, describe:

how the spectral solution differs from the upwind solution,
how the spectral solution differs from the central-difference RK4 solution,
whether the bell-shaped anomaly maintains its original shape,
why cyclic boundary conditions are required for spectral methods.

Comparison Discussion

Compare the results from:

4_3-advection3_central_rk4.png
4_3-advection3_spectral.png

Discuss:

Which scheme best preserves the tracer shape?
Which scheme appears most accurate?
Which scheme would you choose for a smooth periodic tracer field?
What trade-offs exist between simplicity, stability, and accuracy?

Optional Further Experiments

Students interested in further exploration may investigate one or more of the following:

Change the advection velocity \(U\).
Change the grid spacing \(DX\).
Change the time step \(DT\).
Change the initial tracer width \(K\).
Compare Forward Euler and RK4 for the spectral scheme.
Compare all four combinations:
Upwind + Euler
Central + Euler
Central + RK4
Spectral + RK4
Examine how numerical diffusion and numerical instability depend on spatial and temporal resolution.

Discussion Questions

Why does the upwind scheme remain stable even though it is less accurate?
Why can a highly accurate spatial derivative still become unstable?
Why does improving the time integration method improve the solution?
What are the advantages and disadvantages of spectral methods?
Why are cyclic boundary conditions naturally suited for Fourier-based methods?
How do the lessons learned in this chapter apply to weather and climate models?