5.2 Advection-Diffusion Equation

Learning Objectives

After completing this section, students should be able to:

combine advection and diffusion tendencies within a single model
understand the physical roles of transport and mixing
compare advection-diffusion solutions using different advection schemes
evaluate the impact of numerical diffusion and physical diffusion
investigate the effects of time integration schemes on numerical stability
compare finite-difference and spectral advection methods

Concept

Advection and Diffusion Together

Many weather and climate models include multiple physical processes that act simultaneously on the model state.

In this section, we combine advection and diffusion into a single equation:

\[ \frac{\partial T}{\partial t} = -u \frac{\partial T}{\partial x} + \frac{\partial}{\partial x} \left( K \frac{\partial T}{\partial x} \right). \]

The left-hand side represents the time tendency of the tracer.

The right-hand side consists of two physical tendencies:

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

and

\[ F_{\rm diff} = \frac{\partial}{\partial x} \left( K \frac{\partial T}{\partial x} \right). \]

Therefore,

\[ \frac{\partial T}{\partial t} = F_{\rm adv} + F_{\rm diff}. \]

Physical Meaning of Each Term

The first term,

\[ -u \frac{\partial T}{\partial x}, \]

represents advection.

Advection transports the tracer downstream without directly changing its shape. The transport speed is controlled by the velocity \(u\).

The second term,

\[ \frac{\partial}{\partial x} \left( K \frac{\partial T}{\partial x} \right), \]

represents diffusion.

Diffusion acts to reduce spatial gradients by transporting tracer from regions of high concentration toward regions of low concentration. The strength of the mixing is controlled by the diffusion coefficient \(K\).

Physical Interpretation

For a bell-shaped tracer anomaly:

advection moves the anomaly downstream,
diffusion broadens the anomaly,
diffusion reduces the peak amplitude.

Thus,

advection controls transport,
diffusion controls mixing,
together they determine the evolution of the tracer field.

Many weather and climate models include both processes because atmospheric and oceanic tracers are continuously transported by large-scale circulation while simultaneously mixed by smaller-scale turbulent motions.

Numerical Considerations

The advection schemes introduced in Chapter 4 have different numerical properties:

Scheme	Characteristics
Upwind	Stable but diffusive
Central Difference	Less diffusive but can be unstable
Central Difference + RK4	More accurate and more stable
Spectral + RK4	Highly accurate for smooth periodic fields

Because physical diffusion is now included, some numerical oscillations may be damped naturally. However, the choice of advection scheme can still strongly affect the solution.

In this section, we compare several advection schemes while using the same diffusion formulation.

Model / Equation

Governing Equation

The one-dimensional advection-diffusion equation is

\[ \frac{\partial T}{\partial t} = -u \frac{\partial T}{\partial x} + \frac{\partial}{\partial x} \left( K \frac{\partial T}{\partial x} \right). \]

The first term on the right-hand side represents advection,

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

while the second term represents diffusion,

\[ F_{\rm diff} = \frac{\partial}{\partial x} \left( K \frac{\partial T}{\partial x} \right). \]

Therefore,

\[ \frac{\partial T}{\partial t} = F_{\rm adv} + F_{\rm diff}. \]

Discrete System

As in the previous sections, we write the tendency equation as

\[ \frac{\partial T}{\partial t} = F_{\rm adv} + F_{\rm diff}. \]

Using Forward Euler time integration,

\[ T_i^{n+1} = T_i^n + \Delta t \left( F_{\rm adv} + F_{\rm diff} \right)_i. \]

The diffusion tendency is computed using

\[ F_{\rm diff} = K \frac{ T_{i+1} - 2T_i + T_{i-1} } {\Delta x^2}. \]

or equivalently through the flux-divergence form introduced in Section 5.1.

Because diffusion tends to smooth spatial gradients, the diffusion term is generally numerically stable. Therefore, in this section we continue to use the Forward Euler method for the diffusion update.

The advection tendency is computed separately:

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

Unlike diffusion, the numerical behavior of the advection term depends strongly on the spatial discretization and time integration schemes.

Therefore, we will compare several advection methods:

Advection Scheme	Time Integration
Upwind	Forward Euler
Central Difference	Forward Euler
Central Difference	RK4
Spectral	RK4

The diffusion tendency will be treated identically in all experiments, allowing us to focus on how the advection schemes influence the numerical solution.

Numerical Implementation

In the Python code, the model first computes the advection tendency:

Fadv = advection_tendency(T, x, U)

This corresponds to the Forward Euler estimate of the advection tendency.

For higher-order time integration, the RK4 method can be used instead. In this case, the tendency is evaluated four times:

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)

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

After the advection tendency is computed, the model state is updated using

T = T + DT * Fadv

In the advection-diffusion experiments, the same framework can be extended by adding the diffusion tendency:

Fdiff = diffusion_tendency(T, x, Kdiff)

T = T + DT * (Fadv + Fdiff)

The advection tendency may be computed using either Forward Euler or RK4, while the diffusion tendency will continue to use the simpler Forward Euler approach throughout this section.

Sample Code

The template file for this section is

Template: 5_2-diffusion_template.py

This script includes both advection and diffusion tendencies:

Fadv = advection_tendency(T, x, U)
Fdiff = diffusion_tendency(T, x, Kdiff)

T = T + DT * (Fadv + Fdiff)

The default configuration is

U = 0.0
Kdiff = 0.1

so the model initially behaves as a pure diffusion model.

Exercise

Step 1. From Diffusion to Advection-Diffusion

First, create a working copy of the template script:

cp 5_2-diffusion_template.py 5_2-advection_diffusion.py

Step 1.1 Run the Default Diffusion Case

Run the script without modification:

python 5_2-advection_diffusion.py

The default parameters are

U = 0.0
Kdiff = 0.1

This case includes diffusion only.

Questions:

Does the anomaly remain centered near its initial location?
Does the anomaly broaden with time?
Does the peak amplitude decrease?
Why does the anomaly not move downstream?

Step 1.2 Turn Off Diffusion and Activate Advection

Next, set

U = 0.5
Kdiff = 0.0

In the advection_tendency() function, select the upwind scheme:

scheme = "upwind"
#scheme = "central"
#scheme = "spectral"

Use Forward Euler for the advection tendency:

for n in range(NT - 1):
    # -- Euler --
    Fadv = advection_tendency(T, x, U)

    # -- RK4 --
    #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)
    #Fadv = (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0

    Fdiff = diffusion_tendency(T, x, Kdiff)

    T = T + DT * (Fadv + Fdiff)

    T_save[n + 1, :] = T

Run the model again:

python 5_2-advection_diffusion.py

This case includes advection only.

Questions:

Does the anomaly move downstream?
Does the anomaly broaden?
How does this result compare with the pure diffusion case?
Which term controls the downstream transport?
Can you identify any numerical diffusion associated with the upwind scheme?

Step 1.3 Combine Advection and Diffusion

Keep the upwind advection scheme:

scheme = "upwind"
#scheme = "central"
#scheme = "spectral"

Keep Forward Euler for the advection tendency:

for n in range(NT - 1):
    # -- Euler --
    Fadv = advection_tendency(T, x, U)

    # -- RK4 --
    #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)
    #Fadv = (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0

    Fdiff = diffusion_tendency(T, x, Kdiff)

    T = T + DT * (Fadv + Fdiff)

    T_save[n + 1, :] = T

Then activate both processes:

U = 0.5
Kdiff = 0.2

Run the model again:

python 5_2-advection_diffusion.py

This case includes both transport and smoothing.

Questions:

Does the anomaly move downstream?
Does the anomaly broaden with time?
Does the peak amplitude decrease?
How does this solution differ from Step 1.2?
Which changes are caused by advection?
Which changes are caused by diffusion?

Step 2. Central-Difference Advection with Diffusion

In this step, we replace the upwind advection scheme with the central-difference scheme.

Create a new script:

cp 5_2-advection_diffusion.py 5_2-advection_diffusion_central.py

In the advection_tendency() function, select the central-difference scheme:

#scheme = "upwind"
scheme = "central"
#scheme = "spectral"

Continue to use Forward Euler time integration:

for n in range(NT - 1):
    # -- Euler --
    Fadv = advection_tendency(T, x, U)

    # -- RK4 --
    #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)
    #Fadv = (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0

    Fdiff = diffusion_tendency(T, x, Kdiff)

    T = T + DT * (Fadv + Fdiff)

    T_save[n + 1, :] = T

Use

U = 0.5

for all experiments below.

Step 2.1 Central Difference without Diffusion

First, turn off diffusion:

Kdiff = 0.0

Run the model:

python 5_2-advection_diffusion_central.py

Questions:

Does the solution remain stable?
Do oscillations develop near the tracer peak?
How does this result compare with the upwind scheme?
Why is the central-difference scheme more susceptible to numerical instability?

Step 2.2 Central Difference with Strong Diffusion

Next, increase the diffusion coefficient:

Kdiff = 0.2

Continue to use the central-difference scheme and Forward Euler time integration.

Run the model:

python 5_2-advection_diffusion_central.py

Questions:

Does the solution remain stable throughout the simulation?
Are most of the oscillations suppressed?
How does the tracer shape compare with the upwind solution?
Does the larger diffusion coefficient appear to compensate for the instability of the central-difference scheme?

Step 2.3 Central Difference with Weak Diffusion

Now reduce the diffusion coefficient:

Kdiff = 0.1

Continue to use the central-difference scheme and Forward Euler time integration.

Run the model again:

python 5_2-advection_diffusion_central.py

Questions:

Do oscillations become more visible?
Does the solution remain stable?
How does the result compare with the stronger diffusion case?
What trade-off exists between numerical stability and excessive smoothing?

Step 2.4 Central Difference with RK4

Keep

Kdiff = 0.1

and replace the Forward Euler advection tendency with RK4:

for n in range(NT - 1):

    # -- RK4 Advection --
    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)

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

    # -- Diffusion --
    Fdiff = diffusion_tendency(T, x, Kdiff)

    T = T + DT * (Fadv + Fdiff)

    T_save[n + 1, :] = T

Run the model:

python 5_2-advection_diffusion_central.py

Questions:

Does RK4 improve the solution?
Are the oscillations reduced relative to Step 2.3?
How does the solution compare with Step 2.2?
What role does the time integration scheme play in numerical stability?

Step 3. Spectral Advection with RK4

In this step, we replace the central-difference advection scheme with the spectral advection scheme.

Create a new script:

cp 5_2-advection_diffusion_central.py \
   5_2-advection_diffusion_spectral.py

In the advection_tendency() function, select the spectral scheme:

#scheme = "upwind"
#scheme = "central"
scheme = "spectral"

Use

U = 0.5
Kdiff = 0.1

and apply RK4 for the advection tendency:

for n in range(NT - 1):

    # -- RK4 Advection --
    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)

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

    # -- Diffusion --
    Fdiff = diffusion_tendency(T, x, Kdiff)

    T = T + DT * (Fadv + Fdiff)

    T_save[n + 1, :] = T

Run the model:

python 5_2-advection_diffusion_spectral.py

Questions:

Does the solution remain stable throughout the simulation?
How does the tracer shape compare with the central-difference RK4 solution?
Does the spectral method preserve the tracer structure more accurately?
How does diffusion modify the spectral solution?

Step 3.2 Sensitivity to the Diffusion Coefficient

Increase the diffusion coefficient:

Kdiff = 0.2

Run the model again.

Questions:

Does the larger diffusion coefficient noticeably broaden the tracer?
How does the solution compare with the \(Kdiff=0.1\) case?
Does diffusion dominate the solution more strongly?
How does the spectral solution compare with the upwind and central-difference solutions?

Homework

Complete the following three experiments and submit the requested output files.

1. Upwind Advection without Diffusion

Use:

scheme = "upwind"

Forward Euler time integration:

Fadv = advection_tendency(T, x, U)

and

U = 0.5
Kdiff = 0.0

Submit:

5_2-advection_diffusion_upwind.mp4
5_2-advection_diffusion_upwind.png

2. Central Difference with RK4 and Weak Diffusion

Use:

scheme = "central"

RK4 time integration:

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)

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

and

U = 0.5
Kdiff = 0.1

Submit:

5_2-advection_diffusion_central_RK4.mp4
5_2-advection_diffusion_central_RK4.png

3. Spectral Advection with RK4 and Weak Diffusion

Use:

scheme = "spectral"

RK4 time integration:

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)

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

and

U = 0.5
Kdiff = 0.1

Submit:

5_2-advection_diffusion_spectral_RK4.mp4
5_2-advection_diffusion_spectral_RK4.png

Comparison Discussion

Compare the three simulations:

Advection Scheme	Time Integration	Kdiff
Upwind	Euler	0.0
Central Difference	RK4	0.1
Spectral	RK4	0.1

In a short discussion (approximately 1–2 pages), address the following questions:

How does the tracer evolution differ among the three experiments?
Which experiment exhibits the strongest numerical diffusion?
Which experiment best preserves the tracer shape?
How does physical diffusion differ from numerical diffusion?
How does RK4 improve the central-difference solution?
What advantages does the spectral method provide?
Which configuration would you choose for a smooth periodic tracer field, and why?

Optional Further Experiments

Students interested in further exploration may investigate:

different values of \(U\),
different values of \(Kdiff\),
higher spatial resolution (\(\Delta x\)),
smaller time steps (\(\Delta t\)),
central-difference advection with Forward Euler,
spectral advection with no diffusion,
sensitivity of the solution to the diffusion coefficient.

These experiments provide insight into how transport, mixing, numerical stability, and numerical accuracy interact within a numerical model.