Chapter 16.3: Turbulence & Health Metrics

16.3.1 The k-epsilon Turbulence Model

While mixing-length theory (Chapter 16.2) provides a useful first approximation for eddy diffusivity, real street canyon flows require a more sophisticated turbulence model. The standard $k\text{-}\varepsilon$ model solves two transport equations: one for turbulence kinetic energy $k = \tfrac{1}{2}\langle u_i' u_i' \rangle$ and one for its dissipation rate $\varepsilon$.

Turbulence Kinetic Energy Equation

$$\frac{\partial k}{\partial t} + u_j \frac{\partial k}{\partial x_j} = \frac{\partial}{\partial x_j}\!\left(\frac{\nu_t}{\sigma_k}\frac{\partial k}{\partial x_j}\right) + P_k - \varepsilon$$

where the production term is:

$$P_k = \nu_t \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}$$

Dissipation Rate Equation

$$\frac{\partial \varepsilon}{\partial t} + u_j \frac{\partial \varepsilon}{\partial x_j} = \frac{\partial}{\partial x_j}\!\left(\frac{\nu_t}{\sigma_\varepsilon}\frac{\partial \varepsilon}{\partial x_j}\right) + C_{1\varepsilon}\frac{P_k}{k}\varepsilon - C_{2\varepsilon}\frac{\varepsilon^2}{k}$$

Turbulent Viscosity and Diffusivity

The turbulent viscosity is computed from the local values of $k$ and $\varepsilon$:

$$\nu_t = C_\mu \frac{k^2}{\varepsilon}$$

The turbulent diffusivity for pollutant transport is:

$$K = \frac{C_\mu\,k^2}{Sc_t\,\varepsilon}$$

Standard Constants

Constant	Value	Physical Meaning
$C_\mu$	0.09	Links k, epsilon to turbulent viscosity
$C_{1\varepsilon}$	1.44	Production of dissipation
$C_{2\varepsilon}$	1.92	Destruction of dissipation
$\sigma_k$	1.0	Turbulent Prandtl number for k
$\sigma_\varepsilon$	1.3	Turbulent Prandtl number for epsilon
$Sc_t$	0.7	Turbulent Schmidt number

16.3.2 Health Exposure Metrics

Pedestrian Inhaled Dose

The health-relevant quantity is not concentration alone, but the inhaled dose:

$$D = V_E \cdot T_{\text{canyon}} \cdot C_{\text{street}}$$

where $V_E$ is the breathing (ventilation) rate and $T_{\text{canyon}}$ is the time spent in the street canyon. Typical values:

Resting: $V_E \approx 7.5$ L/min
Light exercise (walking): $V_E \approx 20$ L/min
Heavy exercise (cycling): $V_E \approx 40$ L/min

A cyclist inhales 5 times the dose of a sedentary person for the same exposure time and concentration, making active travellers in polluted canyons especially vulnerable.

WHO Air Quality Guidelines (2021)

Pollutant	Annual Mean	24-hour Mean
NO₂	10 μg/m³	25 μg/m³
PM₂.₅	5 μg/m³	15 μg/m³
PM₁₀	15 μg/m³	45 μg/m³

Ventilation Time Scale

The characteristic time for the canyon to flush its polluted air is:

$$\tau_{\text{vent}} = \frac{W}{C_D\,u_H}$$

For a 20 m wide canyon with $C_D = 0.005$ and $u_H = 5$ m/s: $\tau_{\text{vent}} = 20/(0.005 \times 5) = 800$ s ≈ 13 minutes. During calm periods ($u_H < 1$ m/s), ventilation times can exceed an hour.

16.3.3 MFG Extension: Pollution-Aware Routing

In Module 12, the Mean-Field Game (MFG) framework for traffic routing used a congestion cost $F(x, \rho)$ in the Hamilton-Jacobi-Bellman equation. We can replace or augment this with a pollution-aware cost:

$$-\frac{\partial u}{\partial t} + H(x, \nabla u) = F(x, \rho) + \lambda\,C_{\text{OSPM}}(x, \rho)$$

where $C_{\text{OSPM}}(x, \rho)$ is the OSPM concentration at street segment $x$ given traffic density $\rho$, and $\lambda$ is a policy parameter controlling the weight of pollution exposure relative to travel time.

The key coupling: traffic density $\rho$ determines emissions (via HBEFA), which feed into OSPM to compute $C_{\text{OSPM}}$. This modifies the value function $u$, which in turn changes the optimal routing policy, which changes $\rho$. The MFG fixed-point now jointly equilibrates congestion and pollution.

$$\rho \xrightarrow{\text{HBEFA}} S_e \xrightarrow{\text{OSPM}} C_{\text{OSPM}} \xrightarrow{\text{HJB}} u \xrightarrow{\text{FPK}} \rho$$

16.3.4 Python: k-epsilon Solver and Pollution Map

We solve the 1D (vertical) steady-state $k\text{-}\varepsilon$ equations iteratively for a street canyon, then compute the turbulent diffusivity and visualise the results. We also apply OSPM to a synthetic street network.

k-epsilon Canyon Solver & Pollution Map

Python

script.py216 lines

import numpy as np
import matplotlib.pyplot as plt

# ─── k-epsilon 1D canyon solver ───
# Solve steady-state k and epsilon across canyon height
# Simplified 1D: d/dz(nu_t/sigma_k dk/dz) + P_k - eps = 0
#                d/dz(nu_t/sigma_eps deps/dz) + C1e*Pk/k*eps - C2e*eps^2/k = 0

# Constants
C_mu = 0.09
C_1e = 1.44
C_2e = 1.92
sigma_k = 1.0
sigma_eps = 1.3
Sc_t = 0.7
kappa = 0.41

# Canyon parameters
H = 20.0          # building height (m)
Nz = 100          # grid points
dz = H / (Nz - 1)
z = np.linspace(0, H, Nz)

# Rooftop conditions
u_H = 5.0         # m/s
u_star = kappa * u_H / np.log(0.3 * H / (0.1 * H))  # friction velocity

# Mean velocity profile (simplified linear in canyon)
U = u_star * z / H  # linear approximation

# Velocity gradient
dUdz = np.gradient(U, dz)

# Initial conditions
k = np.ones(Nz) * 0.5 * u_star**2
eps = np.ones(Nz) * u_star**3 / (kappa * np.maximum(z, dz))
eps = np.clip(eps, 1e-6, 100)

# Boundary conditions
k_roof = 1.5 * (0.05 * u_H)**2  # TI ~ 5% at roof
eps_roof = C_mu**0.75 * k_roof**1.5 / (0.1 * H)

# Iterative solver
for iteration in range(500):
    nu_t = C_mu * k**2 / np.maximum(eps, 1e-10)
    nu_t = np.clip(nu_t, 1e-6, 100)

# Production
    P_k = nu_t * dUdz**2

# k equation: diffusion + production - dissipation
    k_new = k.copy()
    for i in range(1, Nz - 1):
        diff_k = (nu_t[i] / sigma_k) * (k[i+1] - 2*k[i] + k[i-1]) / dz**2
        k_new[i] = k[i] + 0.01 * (diff_k + P_k[i] - eps[i])

# epsilon equation
    eps_new = eps.copy()
    for i in range(1, Nz - 1):
        diff_eps = (nu_t[i] / sigma_eps) * (eps[i+1] - 2*eps[i] + eps[i-1]) / dz**2
        if k[i] > 1e-10:
            src = C_1e * P_k[i] * eps[i] / k[i] - C_2e * eps[i]**2 / k[i]
        else:
            src = 0
        eps_new[i] = eps[i] + 0.01 * (diff_eps + src)

# Boundary conditions
    k_new[0] = 1e-4       # near-wall
    k_new[-1] = k_roof
    eps_new[0] = C_mu**0.75 * k_new[0]**1.5 / (0.5 * dz)
    eps_new[-1] = eps_roof

# Enforce positivity
    k_new = np.maximum(k_new, 1e-6)
    eps_new = np.maximum(eps_new, 1e-6)

k = k_new
    eps = eps_new

# Compute derived quantities
nu_t_final = C_mu * k**2 / eps
K_diff = nu_t_final / Sc_t  # turbulent diffusivity

# ─── Figure: k-epsilon profiles ───
fig, axes = plt.subplots(1, 4, figsize=(18, 6))
fig.patch.set_facecolor("#0f172a")

quantities = [
    (k, "Turbulence KE k (m²/s²)", "#10b981"),
    (eps, "Dissipation ε (m²/s³)", "#f59e0b"),
    (nu_t_final, "Turb. Viscosity ν_t (m²/s)", "#06b6d4"),
    (K_diff, "Turb. Diffusivity K (m²/s)", "#a78bfa"),
]

for ax, (qty, label, color) in zip(axes, quantities):
    ax.plot(qty, z, color=color, linewidth=2.5)
    ax.set_xlabel(label, color="white", fontsize=10)
    ax.set_ylabel("Height z (m)", color="white", fontsize=11)
    ax.set_facecolor("#0f172a")
    ax.tick_params(colors="gray", labelsize=9)
    ax.grid(True, alpha=0.2, color="#334155")
    for spine in ax.spines.values():
        spine.set_color("#334155")
    ax.axhline(H, color="#ef4444", linestyle="--", alpha=0.5, linewidth=1)

axes[0].set_title("k-ε Profiles Across Canyon", color="white", fontsize=12, fontweight="bold")
plt.tight_layout()
plt.savefig("output.png", dpi=120, bbox_inches="tight", facecolor="#0f172a")
plt.close()

# ─── Figure 2: Grid Street Network Pollution Map ───
fig2, ax2 = plt.subplots(1, 1, figsize=(10, 10))
fig2.patch.set_facecolor("#0f172a")

np.random.seed(42)

# Build a realistic grid street network (8x6 blocks with perturbation)
n_rows, n_cols = 7, 9
block_w, block_h = 120, 90  # meters
jitter = 8  # slight perturbation for realism

# Generate intersection coordinates
nodes = {}
for r in range(n_rows):
    for c in range(n_cols):
        x = c * block_w + np.random.uniform(-jitter, jitter)
        y = r * block_h + np.random.uniform(-jitter, jitter)
        nodes[(r, c)] = (x, y)

# Draw edges and compute OSPM for each street segment
edges = []
# Horizontal streets
for r in range(n_rows):
    for c in range(n_cols - 1):
        edges.append(((r, c), (r, c + 1)))
# Vertical streets
for r in range(n_rows - 1):
    for c in range(n_cols):
        edges.append(((r, c), (r + 1, c)))

for (n1, n2) in edges:
    x0, y0 = nodes[n1]
    x1, y1 = nodes[n2]
    seg_len = np.sqrt((x1 - x0)**2 + (y1 - y0)**2)

# Assign canyon properties (vary by location — denser downtown)
    cx = (x0 + x1) / 2
    cy = (y0 + y1) / 2
    dist_to_centre = np.sqrt((cx - 480)**2 + (cy - 270)**2) / 400
    # Downtown: taller, narrower, more traffic
    H_local = 25 + 30 * max(0, 1 - dist_to_centre) + np.random.uniform(-3, 3)
    W_local = 12 + 10 * min(1, dist_to_centre) + np.random.uniform(-2, 2)
    alpha_local = H_local / max(W_local, 5)
    traffic_local = 1800 * max(0.2, 1 - 0.6 * dist_to_centre) + np.random.uniform(-200, 200)
    traffic_local = max(100, traffic_local)

# Simplified OSPM
    NOx_gkm = 0.08  # Euro 6 diesel
    Q_line = NOx_gkm * 1000 * traffic_local / 3600  # ug/(m*s)
    u_eff = max(0.1, u_H * 0.15 / (1 + alpha_local * 0.3))
    C_D_conc = Q_line / (W_local * u_eff)
    C_R_conc = C_D_conc * min(alpha_local, 3.0) * 0.8
    C_total = C_D_conc + C_R_conc + 15.0  # background

# Color by WHO risk
    if C_total < 25:
        color = "#22c55e"
    elif C_total < 40:
        color = "#eab308"
    elif C_total < 60:
        color = "#f97316"
    else:
        color = "#ef4444"

lw = 1.5 + min(C_total / 60.0, 1.0) * 3
    ax2.plot([x0, x1], [y0, y1], color=color, linewidth=lw, alpha=0.85, solid_capstyle='round')

# Draw intersection nodes
for (r, c), (x, y) in nodes.items():
    ax2.plot(x, y, 'o', color='#94a3b8', markersize=2.5, zorder=5)

# Legend
for label, c in [("< 25 ug/m3 (Good)", "#22c55e"), ("25-40 ug/m3 (Moderate)", "#eab308"),
                  ("40-60 ug/m3 (Unhealthy)", "#f97316"), ("> 60 ug/m3 (Hazardous)", "#ef4444")]:
    ax2.plot([], [], color=c, linewidth=3, label=label)
ax2.legend(fontsize=10, facecolor="#1e293b", edgecolor="#334155", labelcolor="white", loc="upper right")

ax2.set_title("Street Network Pollution Map (OSPM)", color="white", fontsize=14, fontweight="bold")
ax2.set_xlabel("x (m)", color="white")
ax2.set_ylabel("y (m)", color="white")
ax2.set_facecolor("#0f172a")
ax2.tick_params(colors="gray")
ax2.set_aspect("equal")
for spine in ax2.spines.values():
    spine.set_color("#334155")

plt.tight_layout()
plt.savefig("output.png", dpi=120, bbox_inches="tight", facecolor="#0f172a")
plt.close()

# ─── Print health metrics ───
print("=== Health Exposure Analysis ===")
print(f"Canyon: H = {H} m, u_H = {u_H} m/s")
print(f"k-ε converged: k_max = {k.max():.4f} m²/s², K_max = {K_diff.max():.4f} m²/s")
print()
print("─── Ventilation Times ───")
for W_val in [10, 15, 20, 30]:
    for u_val in [1, 3, 5, 8]:
        tau = W_val / (0.005 * u_val)
        print(f"W = {W_val:2d} m, u_H = {u_val} m/s → τ_vent = {tau:6.0f} s ({tau/60:.1f} min)")
print()
print("─── Pedestrian Dose (1 hr exposure at 40 μg/m³) ───")
for activity, V_E in [("Resting", 7.5), ("Walking", 20), ("Cycling", 40)]:
    dose = V_E * 60 * 40e-6  # mg
    print(f"{activity:>10s}: V_E = {V_E:4.1f} L/min → dose = {dose:.3f} mg/hr")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

16.3.5 Key Takeaways

The $k\text{-}\varepsilon$ model provides spatially varying turbulent diffusivity across the canyon, capturing the physics missed by simple mixing length.
Standard constants ($C_\mu = 0.09, C_{1\varepsilon} = 1.44, C_{2\varepsilon} = 1.92$) are well-validated for urban flow applications.
Health exposure depends on dose = breathing rate x time x concentration; cyclists are 5x more exposed than sedentary persons.
Ventilation time $\tau = W/(C_D u_H)$ can exceed an hour during calm conditions, creating persistent pollution episodes.
The MFG framework from Module 12 can incorporate OSPM pollution costs, jointly optimising for congestion and exposure.

← Advection-Diffusion & OSPM Module 17: SUMO Microsimulation →

Constant	Value	Physical Meaning
\(C_\mu\)	0.09	Links k, epsilon to turbulent viscosity
\(C_{1\varepsilon}\)	1.44	Production of dissipation
\(C_{2\varepsilon}\)	1.92	Destruction of dissipation
\(\sigma_k\)	1.0	Turbulent Prandtl number for k
\(\sigma_\varepsilon\)	1.3	Turbulent Prandtl number for epsilon
\(Sc_t\)	0.7	Turbulent Schmidt number