Gaussian Plumes: Simulation

---

Lecture 11

October 5, 2026

Overview

• Review of Last Class
• Gaussian Plume Model Behavior
• Key Takeaways and Upcoming Schedule
• References

Review of Last Class

The Gaussian Plume Model

\[\begin{aligned} C(x,y,z) &= \frac{Q}{2\pi u\sigma_y\sigma_z}\exp\left(\frac{-y^2}{2\sigma_y^2}\right) \times \\ &\left[\exp\left(\frac{-(z-H)^2}{2\sigma_z^2}\right) + \exp\left(\frac{-(z+H)^2}{2\sigma_z^2}\right)\right] \end{aligned}\]

Number of key assumptions: each was used in converting the advection-diffusion equation to a 2D heat equation and obtaining this solution. Without these holding, this formula does not work.

Questions?

Text: VSRIKRISH to 22333

URL: https://pollev.com/vsrikrish

See Results

Gaussian Plume Model Behavior

Visualizing a 3D Model

3D visualization is challenging. We can only ever look at 2D slices of \(C(x, y, z)\), such as:

• Ground-level view, \(z = 0\).
• Vertical section, \(y = 0\).
• Crosswind profile, fixed \(x\).

The Ground-Level Footprint

downwind_positions = 0.05:0.02:8.0
crosswind_positions = -700:5:700
footprint = [plume_ugm3(x, y, 0.0, plume_params) for y in crosswind_positions, x in downwind_positions]

p_foot = contourf(downwind_positions, crosswind_positions, footprint,
    color=:viridis, linewidth=0, colorbar_title="\nC  [µg/m³]", colorbar_titlefontsize=18, clims=(0, 2000),
    xlabel="Distance downwind  [km]", ylabel="Crosswind  [m]")
plot!(p_foot, size=(1050, 340), left_margin=12mm, right_margin=22mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 1: Ground-level concentration. The domain is 8 km long and 1.4 km wide, so the plume is far narrower than it looks here.

Stability class F: moderately stable, a clear night.

Vertical Section

heights = 0:2:200
section_positions = 0.2:0.02:8.0
section = [plume_ugm3(x, 0.0, z, plume_params) for z in heights, x in section_positions]

# Immediately above the stack the plume is a few metres across and the
# concentration there is enormous. Cap the scale at the ground-level range so
# the descent to the ground stays visible instead of being washed out.
p_section = contourf(section_positions, heights, section, color=:viridis, linewidth=0, colorbar_title="\nC  [µg/m³]", colorbar_titlefontsize=18,
    clims=(0, 2000), xlabel="Distance downwind  [km]", ylabel="Height  [m]")

# Height of the highest concentration at each distance, on a finer vertical grid
fine_heights = 0:0.5:200
max_heights = [fine_heights[argmax([plume_ugm3(x, 0.0, z, plume_params) for z in fine_heights])]
    for x in section_positions]
plot!(p_section, section_positions, max_heights, color=cb_vermillion, linewidth=3,
    label="height of maximum")
hline!(p_section, [plume_params.H], color=:black, linewidth=2, linestyle=:dash,
    label="stack height", legend=:topright, xlims=(0.2, 8), ylims=(0, 200))
plot!(p_section, size=(1050, 430), left_margin=12mm, right_margin=22mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 2: Vertical section through the centreline. The plume leaves the stack at 41 m and reaches the ground by roughly 1 km.

Crosswind Profiles

profile_distances = [1, 3, 8]                  # [km]
profile_colours = [cb_vermillion, cb_blue, cb_green]
profile_styles = [:solid, :dash, :dot]
crosswind = -600:2:600                          # [m]

p_cross = plot(xlabel="Ground-Level Crosswind Distance  [m]", ylabel="C  [µg/m³]",
    legend=:topright)
for i in 1:length(profile_distances)
    xk = profile_distances[i]
    plot!(p_cross, crosswind, [plume_ugm3(xk, y, 0.0, plume_params) for y in crosswind],
        color=profile_colours[i], linestyle=profile_styles[i], label="$(xk) km")
end
plot!(p_cross, size=(1200, 450), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 3: Ground-level concentration across the plume at three distances downwind.

Each curve is a slice across the plume at one distance. Close to the stack the plume is narrow and has barely reached the ground; further out it is wider and flatter, because the same amount of pollutant is spread across more air. The width is sigma y, and it changes by a factor of about fifty across this domain. Keep that in mind for later: any grid fine enough to see the plume far away is too coarse to see it near the stack.

Where Is The Worst Pollution?

centerline_positions = 0.05:0.01:8.0
centerline = [plume_ugm3(x, 0.0, 0.0, plume_params) for x in centerline_positions]
peak_value, peak_index = findmax(centerline)
peak_distance = centerline_positions[peak_index]

p_line = plot(centerline_positions, centerline, color=:black, linewidth=4,
    xlabel="Distance downwind  [km]", ylabel="C  [µg/m³]",
    legend=false)
scatter!(p_line, [peak_distance], [peak_value], color=cb_vermillion, markersize=9)
plot!(p_line, size=(1050, 330), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 4: Ground-level centreline concentration against downwind distance.

Why Does This Occur Downwind?

\[C(x, 0, 0) = \frac{Q}{\pi u\,\sigma_y\sigma_z}\exp\left(-\frac{H^2}{2\sigma_z^2}\right)\]

Near the stack, which factor controls \(C\)? Far downwind?

\(C\) peaks where arrival (the exponential) rises as fast as dilution falls, in percentage terms: \(\sigma_z \approx H/\sqrt{2} =\) 29 m, about 3.3 km, where the exponential is \(e^{-1}\). The actual peak is at 2.5 km: \(\sigma_y\) grows faster than \(\sigma_z\).

Board item. Near the stack sigma z is small compared with H, so the exponential is tiny: the plume has not yet spread down to the ground. Far downwind the exponential is close to 1, and the one over sigma y sigma z factor takes over: the plume is diluted into more and more air. Call these arrival, the exponential, and dilution, the one over sigma y sigma z factor. The peak sits where the two effects balance in rate, not in size. Moving downwind, each 1% growth in sigma z cuts the dilution factor by 2%, if sigma y grows in step with sigma z, and raises the exponential by H squared over sigma z squared percent, which is large near the stack and shrinks with distance. They match when H squared over sigma z squared is 2, which is exactly where the exponential equals e to the minus 1. So the e to the minus 1 comes from the 2 in the dilution factor, not from anything special about the exponential. Derivation for the board: assume sigma y is a fixed multiple of sigma z. Then C is proportional to one over sigma z squared, times the exponential. Take the log, differentiate with respect to sigma z and set it to zero: minus two over sigma z plus H squared over sigma z cubed equals zero, so sigma z squared equals H squared over two. Here sigma y grows faster than sigma z: at the peak, a 1% increase in sigma z comes with a 1.69% increase in sigma y, so dilution falls by 2.69% and the exponential has to be climbing that fast too, which happens closer to the stack. The same argument predicts the exponent at the peak is half of that, 1.34, and the model’s exponent there is 1.34: the exponential is at about 0.26 rather than 0.37.

Stack Height Changes The Ground-Level \(C\)

standard = 196                                  # 1-hour SO₂ standard, 75 ppb  [µg/m³]
stack_heights = [41, 60, 100]                   # effective stack height  [m]
height_colours = [cb_vermillion, cb_blue, cb_green]
height_styles = [:solid, :dash, :dot]
far_positions = 0.05:0.01:20                    # [km]

height_peaks = zeros(length(stack_heights))
height_peak_distances = zeros(length(stack_heights))
p_height = plot(xlabel="Distance downwind  [km]", ylabel="Centreline C  [µg/m³]",
    legend=:topright)
for i in 1:length(stack_heights)
    params_h = (; plume_params..., H=stack_heights[i])
    centreline_h = [plume_ugm3(x, 0.0, 0.0, params_h) for x in far_positions]
    height_peaks[i], j = findmax(centreline_h)
    height_peak_distances[i] = far_positions[j]
    plot!(p_height, far_positions, centreline_h, color=height_colours[i],
        linestyle=height_styles[i], label="H = $(stack_heights[i]) m")
end
hline!(p_height, [standard], color=:black, linestyle=:dashdot, linewidth=2, label="1-hour SO₂ standard, $(standard) µg/m³")
plot!(p_height, size=(1200, 450), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 5: Ground-level centreline concentration for three effective stack heights, with the same emissions. Dash-dot: the 1-hour SO₂ standard.

Stack height enters through the exponential, so it has an enormous effect on the ground-level peak: going from 41 to 100 metres cuts it by a factor of about sixteen. But the emissions are unchanged. The pollution has not gone away; it comes down further downwind, more dilute, over a larger area. That is the tall-stacks trade-off, and part of how acid rain became a regional problem rather than a local one. H here is the effective height, which includes how far the hot exhaust rises above the physical stack. One caution: the 100 m peak is at nearly 15 km, and the fitted sigma curves are not reliable much beyond 10 km, so treat that number as rough.

Stability Changes The Profile Too

# Martin (1976) coefficients: sigma_y = a x^0.894, sigma_z = c x^d + f, with x in km
stability_coefs = Dict(
    "B" => (; a=156, c_near=106.6, d_near=1.149, f_near=3.3, c_far=108.2, d_far=1.098, f_far=2),
    "D" => (; a=68, c_near=33.2, d_near=0.725, f_near=-1.7, c_far=44.5, d_far=0.516, f_far=-13),
    "F" => (; a=34, c_near=14.35, d_near=0.740, f_near=-0.35, c_far=62.6, d_far=0.180, f_far=-48.6))

# Ground-level centreline concentration [µg/m³] under a given stability class
function centreline_class(xk, coef, params)
    if xk <= 0.05
        return 0.0
    end
    sy = coef.a * xk^0.894
    if xk < 1
        sz = coef.c_near * xk^coef.d_near + coef.f_near
    else
        sz = coef.c_far * xk^coef.d_far + coef.f_far
    end
    return params.Q / (π * params.u_wind * sy * sz) * exp(-params.H^2 / (2sz^2)) * 1e6
end

class_names = ["B", "D", "F"]
class_labels = ["B (sunny afternoon)", "D (overcast)", "F (clear night)"]
class_colours = [cb_orange, cb_purple, cb_blue]
class_styles = [:solid, :dash, :dot]
near_positions = 0.05:0.005:10                  # [km]

class_peaks = zeros(length(class_names))
class_peak_distances = zeros(length(class_names))
p_class = plot(xlabel="Distance downwind  [km]", ylabel="C  [µg/m³]", legend=:topright)
for i in 1:length(class_names)
    c = [centreline_class(x, stability_coefs[class_names[i]], plume_params) for x in near_positions]
    class_peaks[i], j = findmax(c)
    class_peak_distances[i] = near_positions[j]
    plot!(p_class, near_positions, c, color=class_colours[i], linestyle=class_styles[i],
        label=class_labels[i])
end
hline!(p_class, [standard], color=:black, linestyle=:dashdot, linewidth=2, label="1-hour SO₂ standard, $(standard) µg/m³")
plot!(p_class, size=(1200, 450), left_margin=12mm, right_margin=8mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 6: Ground-level centreline concentration for the same stack under three stability classes. Dash-dot: the 1-hour SO₂ standard.

Same stack, same emissions, same wind; only the stability changes. On a sunny afternoon the ground heats the air and strong vertical eddies drag the plume down within a few hundred metres, so the peak is high and close to the stack. On a clear night the air is stable, the plume stays aloft as a thin ribbon and comes down slowly, so the peak is lower but much further out, and the plume stays above the standard for many kilometres. So which class is the worst case depends on the question: the highest peak, or the most ground above the standard, or a particular receptor. You will work through that on HW5. Holding the wind at 2.5 m/s for all three is a simplification, but that speed is plausible for both a sunny afternoon and a clear night.

The Area Above the Standard

area_dx, area_dy = 0.02, 2                      # grid spacing: [km] downwind, [m] crosswind
area_x = (area_dx/2):area_dx:10
area_y = (-1000 + area_dy/2):area_dy:1000
area_C = [plume_ugm3(x, y, 0.0, plume_params) for y in area_y, x in area_x]
# Each grid cell above the standard contributes its own area, dx × dy.
area_above = count(area_C .> standard) * (area_dx * 1000) * area_dy / 1e6   # [km²]

p_area = contourf(area_x, area_y, area_C, color=:viridis, linewidth=0, colorbar_title="\nC  [µg/m³]", colorbar_titlefontsize=18, clims=(0, 2000),
    xlabel="Distance downwind  [km]", ylabel="Crosswind  [m]")
contour!(p_area, area_x, area_y, area_C, levels=[standard], color=:black, linewidth=3,
    colorbar_entry=false)
plot!(p_area, size=(1050, 310), left_margin=12mm, right_margin=22mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 7: Ground-level footprint with the standard’s contour (black). The exceedance continues past the 10 km edge of the domain.

About 5.0 km² is above the 1-hour SO\(_2\) standard (196 µg/m³) within 10 km; the plume is still above it at the edge.

This is the third question a permit asks: not just how bad the worst point is, but how much ground is exposed above the limit. We count the grid cells above the standard and add up their areas. On a stable night like this one the plume stays thin and travels a long way, so the exceedance runs right off the edge of our domain. The honest answer to “how large is the area” is “at least five square kilometres, and we cannot say how much more from these sigma curves”, because they are not reliable much beyond 10 km.

Impact of the Spatial Grid

# Crosswind samples offset by half a cell, so none sits exactly on the centreline.
# The plume is symmetric, so sample one side and double the area.
function grid_summary(crosswind_spacing; downwind_spacing=0.05)
    xs = (downwind_spacing/2):downwind_spacing:10
    ys = (crosswind_spacing/2):crosswind_spacing:1000
    samples = [plume_ugm3(x, y, 0.0, plume_params) for y in ys, x in xs]
    area = 2 * count(samples .> standard) * (downwind_spacing * 1000) * crosswind_spacing / 1e6
    return maximum(samples), area
end

println("| Crosswind spacing (m) | Peak (µg/m³) | Area above standard (km²) |")
println("|---:|---:|---:|")
for spacing in [400, 100, 20, 2]
    grid_peak, grid_area = grid_summary(spacing)
    println("| $(spacing) | $(Int(round(grid_peak))) | $(round(grid_area; digits=2)) |")
end

Crosswind spacing (m)	Peak (µg/m³)	Area above standard (km²)
400	571	5.4
100	1454	5.05
20	1732	5.03
2	1747	5.03

Only the peak needs a sample near the centreline.

The samples themselves are exact; what changes with the grid is what we compute from them. The peak is a single point, so a coarse grid simply misses it: at 400 metres no sample comes close enough to the centreline, and we report about a third of the true peak. The area is an integral, and a coarse grid’s errors at the edges of the footprint partly cancel, so it is within a few percent even at 400 metres. This is the same recipe as Week 4: decide what number you need to report, then halve the spacing until that number stops moving. Refine for the quantity you report, not for everything.

How To Reduce Pollution?

The 1-hour SO\(_2\) standard is 75 ppb, about 196 µg/m³. At \(H = 41\) m the peak is 1747 µg/m³: about 9 times the standard.

Two ways to comply:

• Raise the stack: 100 m is enough here.
• Cut emissions by a factor of about 9, since \(C\) is proportional to \(Q\).

Some caveats to say out loud. The standard is not a single-hour cap: a site complies if the 99th percentile of each day’s highest one-hour value, averaged over three years, is at or below 75 ppb. Our sigma curves correspond to roughly ten-minute to one-hour averages, so the comparison is about right in averaging time but not in form. And a real permit analysis adds the background concentration already in the air, and models every hour of several years of weather. So this says the source would be a problem on a stable night with this wind, not that it violates the standard as defined.

Impact of Multiple Sources

stacks = [(0, -150), (1, 0), (2, 150)]          # (downwind [km], crosswind [m]) of each stack
stack_params = (; plume_params..., Q=10)        # our stack, with emissions cut tenfold

# The plume equation is linear in Q, so concentrations from separate stacks add.
combined(x, y) = sum(plume_ugm3(x - s[1], y - s[2], 0.0, stack_params) for s in stacks)

super_x = -0.5:0.02:12
super_y = -700:5:700
super_C = [combined(x, y) for y in super_y, x in super_x]
combined_peak, k = findmax(super_C)
combined_peak_distance = super_x[k[2]]
single_peak = maximum(plume_ugm3(x, 0.0, 0.0, stack_params) for x in super_x)

p_super = contourf(super_x, super_y, super_C, color=:viridis, linewidth=0, colorbar_title="\nC  [µg/m³]", colorbar_titlefontsize=18,
    xlabel="Distance downwind  [km]", ylabel="Crosswind  [m]")
contour!(p_super, super_x, super_y, super_C, levels=[standard], color=:black, linewidth=3,
    colorbar_entry=false)
scatter!(p_super, [s[1] for s in stacks], [s[2] for s in stacks], color=:white,
    markerstrokecolor=:black, markersize=9, label=false)
plot!(p_super, size=(1050, 330), left_margin=12mm, right_margin=22mm,
    bottom_margin=12mm, top_margin=5mm)

Figure 8: Ground-level footprint from three facilities like ours (white markers), each cut to 10 g/s. Black: the standard’s contour.

Each facility cuts to 10 g/s and complies alone (175 µg/m³). Together: 282 µg/m³ at 5.5 km, where no single plume peaks.

This picks up the emissions-cutting option from the previous slides. Suppose three facilities like ours (41 m, class F, 2.5 m/s) each cut emissions tenfold, to 10 g/s, so each one complies when checked on its own. The point is that compliance source by source does not add up to compliance overall. The plume equation is linear in the emission rate, so the concentration from several sources is just the sum of their individual plumes. Each of these three stacks would pass a single-source check. But downwind their plumes have spread wide enough to overlap, and the combined maximum lands further downwind than any one stack’s own peak, between their centrelines. This is why permit modelling includes every nearby source plus the background, not each stack on its own.

Key Takeaways and Upcoming Schedule

Key Takeaways

• The ground-level peak is downwind of the stack: the plume must spread down to the ground first.
• A taller stack lowers the peak and moves it out.
• Whether the grid resolution matters depends on what you report.
• Concentrations add and can change the area as well as peak level.

Next Classes

Wednesday: Model calibration and validation.

Assessments

Mini-Project 1: Due Thursday, October 22.

Homework 4: Due Thursday, October 8.

References