5.1 Diffusion Equation

Learning Objectives

understand the physical meaning of diffusion
distinguish between advection and diffusion
implement a diffusion tendency in a numerical model
understand the role of the diffusion coefficient
investigate the numerical stability of diffusion
compare advection and diffusion processes

Concept

What Is Diffusion?

Diffusion is a physical process that acts to reduce spatial gradients.

When the tracer concentration is larger in one location than nearby regions, diffusion transports material from regions of high concentration toward regions of low concentration.

As a result, diffusion smooths spatial gradients and spreads anomalies over time.

Examples include:

heat conduction
molecular diffusion
turbulent mixing
atmospheric eddy diffusion
ocean mixing

Unlike advection, diffusion does not transport a tracer as a coherent structure. Instead, diffusion broadens and weakens anomalies.

Advection Versus Diffusion

The advection equation is

\[ \frac{\partial T}{\partial t} = -u \frac{\partial T}{\partial x}. \]

Advection transports a tracer.

The diffusion equation is

\[ \frac{\partial T}{\partial t} = K \frac{\partial^2 T}{\partial x^2}. \]

Diffusion smooths a tracer.

For a bell-shaped anomaly:

advection shifts the anomaly downstream
diffusion broadens the anomaly
diffusion reduces the peak amplitude

Model / Equation

Governing Equation

In this chapter, we solve the one-dimensional diffusion equation:

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

where the diffusion tendency is defined as

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

Here,

\(T\) is the tracer,
\(K\) is the diffusion coefficient,
\(F_{\rm diff}\) is the diffusion tendency,
\(x\) is the spatial coordinate,
\(t\) is time.

This form shows that diffusion is related to the divergence of a diffusive flux. The inner derivative,

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

represents the spatial gradient of the tracer. The diffusion coefficient \(K\) controls the strength of the mixing. The outer derivative,

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

describes how the diffusive flux changes in space.

For constant \(K\), the equation reduces to

\[ F_{\rm diff} = K \frac{\partial^2 T}{\partial x^2}. \]

Thus, the one-dimensional diffusion equation can also be written as

\[ \frac{\partial T}{\partial t} = K \frac{\partial^2 T}{\partial x^2}. \]

This equation indicates that diffusion depends on the curvature of the tracer field. Regions with strong curvature experience stronger diffusion tendencies.

Discrete System

As in the advection chapter, the diffusion equation is written as

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

Using Forward Euler time integration,

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

The diffusion tendency is

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

For simplicity, we first assume that \(K\) is constant.

Define the diffusive flux as

\[ G = K \frac{\partial T}{\partial x}. \]

Using a finite difference,

\[ G_i = K \frac{T_{i+1}-T_i}{\Delta x}, \]

and

\[ G_{i-1} = K \frac{T_i-T_{i-1}}{\Delta x}. \]

The diffusion tendency is obtained from the spatial derivative of this flux:

\[ \left(F_{\rm diff}\right)_i = \frac{G_i-G_{i-1}} {\Delta x}. \]

Substituting the expressions for \(G_i\) and \(G_{i-1}\),

\[ \left(F_{\rm diff}\right)_i = \frac{ K\frac{T_{i+1}-T_i}{\Delta x} - K\frac{T_i-T_{i-1}}{\Delta x} } {\Delta x}. \]

Since \(K\) is constant, this simplifies to

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

Therefore, the Forward Euler update becomes

\[ T_i^{n+1} = T_i^n + \Delta t K \frac{ T_{i+1} - 2T_i + T_{i-1} } {\Delta x^2}. \]

In the exercises, we will implement diffusion in two equivalent ways:

compute the diffusive flux \(G = K\partial T/\partial x\) first, then compute \(\partial G/\partial x\);
use the simplified second-derivative form directly.

When \(K\) is constant, both approaches produce the same numerical solution.

Sample Code

This section builds upon the summarized advection model developed in the previous chapter:

Template: 4_3-advection3_summary.py

The script already includes:

a one-dimensional grid,
a bell-shaped initial condition,
cyclic boundary conditions,
selectable advection schemes,
Forward Euler and RK4 time integration.

In this chapter, we will modify the advection model to include diffusion and eventually combine both processes within a single model.

Exercise

Step 1. Prepare the Diffusion Model

Create a new script from the advection summary model:

cp 4_3-advection3_summary.py 5_1-diffusion.py

Set the advection velocity to zero:

U = 0.0

Then define a diffusion coefficient near the parameter section:

Kdiff = 0.2

Next, create a new diffusion tendency function:

def diffusion_tendency(T, x, Kdiff):

    Fdiff = np.zeros_like(T)

    return Fdiff

In the time integration loop, replace the advection tendency with the diffusion tendency.

For Forward Euler:

Fdiff = diffusion_tendency(T, x, Kdiff)

T = T + DT * Fdiff

Run the model:

python 5_1-diffusion.py

Questions:

Does the bell-shaped anomaly change with time?
Why does nothing happen when Fdiff = 0?
Which part of the code controls the diffusion tendency?
Why did we set U = 0.0 for this first experiment?

This step confirms that the diffusion model framework runs correctly before the diffusion term is implemented.

Step 2. Implement the Diffusion Tendency

In this step, implement the diffusion tendency,

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

For constant \(K\), this becomes

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

Step 2.1 Implement Diffusion Using a Loop

First, implement the diffusion tendency using an explicit loop.

Replace the diffusion_tendency() function with:

def diffusion_tendency(T, x, Kdiff):

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

    for i in range(1, NX-1):
        Fdiff[i] = Kdiff * ( T[i+1] - 2.0 * T[i] + T[i-1]) / DX**2

    # Cyclic boundary points
    Fdiff[0] = Kdiff * ( T[1] - 2.0 * T[0] + T[-1]) / DX**2
    Fdiff[-1] = Kdiff * ( T[0] - 2.0 * T[-1] + T[-2]) / DX**2

    return Fdiff

Run the model:

python 5_1-diffusion.py

Questions:

Does the bell-shaped anomaly move?
Does the anomaly become broader?
Does the peak amplitude decrease?
Why does diffusion smooth the tracer field?

Step 2.2 Implement Diffusion Using `np.diff`

The loop version is useful for understanding the equation, but NumPy can compute differences more efficiently.

Create a new script from the advection summary model:

cp 5_1-diffusion.py 5_1-diffusion_step2.2.py

Replace the diffusion_tendency() function with:

def diffusion_tendency(T, x, Kdiff):

    Fdiff = np.zeros_like(T)

    # First gradient
    G = np.zeros_like(T)
    G[:-1] = Kdiff * np.diff(T) / DX

    # Cyclic boundary for the gradient
    G[-1] = Kdiff * (T[0] - T[-1]) / DX

    # Divergence of the gradient
    Fdiff[1:] = (G[1:] - G[:-1]) / DX

    # Cyclic boundary for the divergence
    Fdiff[0] = (G[0] - G[-1]) / DX

    return Fdiff

This implementation follows the flux-divergence form:

\[ G = K \frac{\partial T}{\partial x}, \]

\[ F_{\rm diff} = \frac{\partial G}{\partial x}. \]

Run the model again.

Questions:

Does this solution match the loop version?
Why are cyclic boundary values needed?
Why is this form useful for more general diffusion problems?
How would this approach help if \(K\) varied in space?

Step 3. Sensitivity to the Diffusion Coefficient

Next, examine how the diffusion coefficient controls the strength of smoothing.

Test the following values:

Kdiff = 0.05

Kdiff = 0.2

Kdiff = 0.5

Run the model for each case.

Questions:

How does increasing Kdiff affect the rate of smoothing?
How does the peak amplitude change?
How does the width of the anomaly change?
Does the total tracer amount appear to be conserved?
What happens if Kdiff is too large?

The diffusion coefficient controls how rapidly gradients are smoothed. Larger values of Kdiff produce stronger mixing.

Key Points

Diffusion acts to reduce spatial gradients and smooth tracer distributions.
Unlike advection, diffusion does not transport a tracer as a coherent structure. Instead, diffusion broadens anomalies and reduces their peak amplitude.
The one-dimensional diffusion equation can be written as

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

The quantity

\[ G = K \frac{\partial T}{\partial x} \]

represents the diffusive flux.

The diffusion tendency is obtained from the divergence of the diffusive flux:

\[ F_{\rm diff} = \frac{\partial G}{\partial x}. \]

For a constant diffusion coefficient \(K\), the diffusion tendency simplifies to

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

Diffusion can be implemented in two equivalent ways:
compute the flux \(G\) first and then compute its divergence,
use the simplified second-derivative form directly.
The flux-divergence form is often preferred because it can easily accommodate spatially varying diffusion coefficients.
Increasing the diffusion coefficient increases the rate at which gradients are smoothed.
Larger diffusion coefficients produce broader tracer distributions and smaller peak amplitudes.
Numerical stability depends on the diffusion number

\[ D = \frac{K \Delta t} {\Delta x^2}. \]

Increasing \(K\), increasing \(\Delta t\), or decreasing \(\Delta x\) makes the numerical solution more susceptible to instability.
Diffusion is commonly used in weather and climate models to represent unresolved mixing processes and to suppress grid-scale numerical noise.

Homework

1. Weak Diffusion Experiment

Use:

Kdiff = 0.05

Run the model and submit:

5_1-diffusion_005.mp4
5_1-diffusion_005.png

Provide a short description addressing:

How does the tracer evolve with time?
How much does the peak amplitude decrease?
How much does the tracer broaden?

2. Strong Diffusion Experiment

Use:

Kdiff = 0.50

Run the model and submit:

5_1-diffusion_050.mp4
5_1-diffusion_050.png

Provide a short description addressing:

How does the tracer evolve with time?
How much does the peak amplitude decrease?
How much does the tracer broaden?

Comparison Discussion

Compare the two experiments:

Kdiff = 0.05
Kdiff = 0.50

Discuss:

Which experiment produces stronger smoothing?
Which experiment reduces the peak amplitude more rapidly?
Which experiment produces a broader tracer distribution?
How does the diffusion coefficient affect the rate of mixing?

Optional Further Experiments

Students interested in additional exploration may investigate:

values of \(Kdiff\) between 0.05 and 0.70,
the effect of changing \(\Delta x\),
the effect of changing \(\Delta t\),
the numerical stability criterion

\[ \frac{K \Delta t} {\Delta x^2} < \frac12, \]

comparison between the flux-divergence form and the second-derivative form of the diffusion tendency.

Discussion Questions

How does diffusion differ from advection?
Why does diffusion reduce gradients?
Why does increasing \(Kdiff\) increase the rate of smoothing?
Why is diffusion commonly included in numerical weather and climate models?
What physical processes might \(Kdiff\) represent in the atmosphere or ocean?

5.1 Diffusion Equation

Learning Objectives

Concept

What Is Diffusion?

Advection Versus Diffusion

Model / Equation

Governing Equation

Discrete System

Sample Code

Exercise

Step 1. Prepare the Diffusion Model

Step 2. Implement the Diffusion Tendency

Step 2.1 Implement Diffusion Using a Loop

Step 2.2 Implement Diffusion Using np.diff

Step 3. Sensitivity to the Diffusion Coefficient

Key Points

Homework

1. Weak Diffusion Experiment

2. Strong Diffusion Experiment

Comparison Discussion

Optional Further Experiments

Discussion Questions

Step 2.2 Implement Diffusion Using `np.diff`