How Plume Calc computes risk

01Model pipeline

1. Inputs. Equipment position, release height, equivalent diameter, efflux velocity and temperature; an operating scenario; aircraft; and an hourly ambient profile of potential temperature, wind speed and, when available, wind direction.
2. Mean-flow plume. Each hour, each source (or merged cluster) is integrated as a turbulent buoyant jet in a crossflow.
3. Aircraft response. For each aircraft and hazard, the vertical gust U_req that causes severe turbulence or upset is computed once.
4. Turbulence model. The mean vertical velocity w at a point defines a gust distribution; the chance of exceeding U_req is the hourly probability there.
5. Aggregation. Hourly probabilities are averaged over valid hours, combined across independent clusters and contoured.

02Mean-flow plume

The plume is a top-hat integral model in its own vertical plane. With arc length s along the centreline, local radius b, speed V and inclination φ above horizontal, the model conserves volume, horizontal momentum, vertical momentum (with buoyancy) and buoyancy flux in a stratified ambient. Per π, the state is

with M = √(M_x² + M_z²), V = M/Q, b = Q/√M, Δθ = F/Q, and governing equations

U(z) is the ambient wind speed, θ_a(z) the ambient potential temperature, dθ_a/dz its gradient (the segment slope of the piecewise-linear profile) and g = 9.80665 m/s². Entrainment is the sum of a shear (jet) term and a crossflow term:

α = 0.11 is the top-hat jet coefficient (Morton, Taylor & Turner 1956 give about 0.08–0.12; Hoult, Fay & Forney 1969 use the same form). β = 0.6 is the crossflow coefficient of Hoult, Fay & Forney (1969).

Release. A vertical discharge (φ = 90°) with b₀ = D/2, V₀ the efflux velocity and Δθ₀ the initial potential-temperature excess. Absolute efflux temperatures are converted to potential temperature using standard-atmosphere pressure at ground elevation plus release height, then reduced by the ambient θ at release height. Equipment given as a temperature excess uses it directly.

Numerics. Fourth-order Runge–Kutta with step ds = clamp(0.2b, 0.05 m, 20 m). Integration stops when the centreline vertical velocity falls below the smallest hazard-relevant velocity (never below 0.05 m/s), the plume passes 3,000 m AGL, the path exceeds 15 km, or 4,000 steps are taken. Ambient θ and U are linearly interpolated on the height grid and held constant outside it.

Gaussian conversion. The turbulence model needs a smooth velocity field, so each top-hat state becomes a Gaussian of equal volume and momentum flux:

Check: πb_g²w_c = πb²V, the top-hat volume flux. w_c is the centreline vertical velocity.

03Severe turbulence

Severe turbulence is an incremental vertical load factor Δn = 1.0 g at approach speed. The required gust follows the FAR 23.341 gust-load formula with the Pratt–Walker alleviation factor:

W/S is weight per wing area (N/m²), c̄ the mean chord, a the lift-curve slope (per radian), V_a the approach speed and ρ₀ sea-level density, so U_req is an equivalent-airspeed gust. The gust is a 1-cosine of length 12.5 chords hitting the aircraft centre.

04Aircraft upset

Upset is a bank-angle change from 25° to 45°, that is 20°. The aircraft is a single-degree-of-freedom roll system with no pilot or autopilot input (after Houbolt & Sen 1971, as used by MITRE). The gust is a 1-cosine of length 12.5 chords acting on one wing only, the worst case of a plume at a wing tip.

S is wing area and B the span. The equation is integrated with RK4 (step the smaller of T/400 and τ/20, where τ = I_xx/L_p) out to T + 10τ, recording the peak bank change. U_req is the gust whose peak equals 20°, found by doubling from 1 m/s to bracket it and then 60 bisection steps. Aircraft with no known I_xx cannot be evaluated for upset and are reported as unavailable, as with MITRE’s Light-Sport reference type.

05Turbulence and probability

Following Papanicolaou & List (1988) as used by MITRE, the RMS gust at a point is proportional to the local mean vertical velocity:

erfc uses the Numerical Recipes Chebyshev fit (fractional error below 1.2×10⁻⁷, which keeps relative accuracy in the tail). Probabilities below 10⁻⁹ (z = 6) are set to zero, so only points with w ≥ U_req / (1 + 0.25·6) = U_req/2.5 matter. That bound limits both the plume integration and the area each trajectory step covers on the grid, a radius of b_g√ln(w_c/w_rel). Where several steps of one trajectory reach a cell in the same hour, the largest local w is used.

06Aggregation

• Side profile. An x–z grid (x downwind, z above ground), auto-sized from the 400 calmest and 150 most unstable hours plus an even spread of about 600 hours (stride coprime with 24), capped at 2,500 m high and 6 km downwind, to roughly 220 × 180 cells with a 0.25 m minimum cell. Each valid hour adds its probability; the result is the mean over valid hours.
• Contours at 10⁻⁷, 10⁻⁶, 10⁻⁵, 10⁻⁴, 10⁻³ and 10⁻², coloured yellow to dark red as in MITRE’s legend.
• Data-sufficiency warning. MITRE recommends three years (about 26,300 hours). A run with fewer than 25,000 valid hours (about 95 %) is flagged, because NOAA data have dropouts. Hours with missing data are excluded, not zero-filled.

07Multiple sources

Clustering. Two sources are linked when both their horizontal separation and their release-height difference are at most k · max(D_a, D_b). Links are transitive (union–find), and k is the merge factor.

Flux-conserving merge. A cluster becomes one equivalent stack with the summed volume and momentum flux. With q_i = π D_i²V_i/4:

This is an AERMOD-style merged stack. Buoyancy is merged each hour because absolute efflux temperatures depend on the ambient.

Combining clusters. Clusters are integrated separately and treated as independent. For each voxel and hour,

and the final value is the mean over hours. Clusters with identical release signatures share one side-profile computation; by default up to 12 profiles are produced, strongest first (ranked by D²V²).

08MITRE-style averaged mode

This reproduces what a single-input tool forces on a multi-source site. All active units become N identical stacks, with mean diameter, velocity, height and absolute temperature (temperature excesses are referenced to a 15 °C day), placed on a line through the site centroid at the mean nearest-neighbour distance of the real layout. They merge only if that averaged spacing is within k diameters. When equipment sits acres apart, averaging changes both what merges and where the plume is placed, so risk can vanish or distort. Use the mode as a comparison, not as the answer.

093D site grid and the risk ceiling

For the map, a voxel grid in local east–north–up metres spans the cluster bounding box plus a margin of max(downwind reach, 60 m) + 30 m. Horizontal cells are at least 5 m and the grid is capped at 2,000,000 voxels. Vertically it has 8 to 60 layers (nominally 7.62 m, or 25 ft) up to 1.3 times the tallest plume reach, with a minimum of 30 m.

Each hour, the side-profile trajectory is swung into plan along the direction the wind blows toward. That direction comes from the wind direction (interpolated on unit vectors) at the mid-height of the trajectory, and the Gaussian is applied radially around the path in 3D. Hours whose plume leaves the grid are counted and reported as clipped. Maps require wind direction, so a MITRE CSV, which has none, yields side profiles only.

The risk ceiling is the highest altitude above each location at which the mean probability still reaches the selected contour level, that is, how high the hazard extends over each point.

10Weather

• Source. NOAA HRRR (ncep_hrrr_conus) hourly data from the Open-Meteo Historical Forecast API, with a global GFS fallback. Requests are chunked into three-month pieces, cached, and retried with exponential back-off.
• Variables. 2 m temperature, 10 m and 80 m wind speed and direction, surface pressure, and temperature, wind speed, wind direction and geopotential height at 1000, 975, 950, 925, 900, 850, 800 and 700 hPa.
• Heights. Level height AGL = geopotential height (MSL) − grid-cell elevation. Levels at or below ground are dropped.
• Potential temperature. θ = T (1000/p)^0.2857, with T in K and p in hPa.
• Resampling. Profiles are linearly resampled to a fixed 40-level grid from 0 to 3,000 m AGL (10–25 m spacing near the ground, 50 m up to 900 m, 100 m up to 2,000 m, 200 m above). Wind direction is interpolated on unit vectors so 350° to 10° passes through 0°. End values are held outside the source range.
• MITRE CSV, both ways. Import and export use MITRE’s format: one row per hour and level with Unix timestamp (s), elevation (m AGL), potential temperature (K) and wind speed (m/s), and a [NO DATA] row for missing hours.

11Units and coordinates

Computation is SI internally (m, s, K, kg); feet and knots are converted at the boundary. Heights are above ground level, and ground elevation (MSL) is used only to set pressure. The map uses WGS84. DXF site plans are read in their own coordinate system, with automatic detection of NAD83 State Plane zones in US survey feet or international feet. NAD83 and WGS84 are treated as identical, which holds to about a metre in the contiguous US.

12Validation

What is tested today:

• Similarity exponents. The plume integrator reproduces the similarity exponents of a pure jet and a pure plume in a uniform, calm ambient.
• Briggs two-thirds law. Buoyant rise in a crossflow follows Briggs’ 2/3 law within ±35 %.
• Gust-loads anchor. The severe-turbulence gust for MITRE’s Light GA example comes out at about 36 ft/s.
• Tail probability. The Q-function is checked against reference values deep into the tail.

Pending. Parity against MITRE tool outputs on identical inputs requires client data and has not been done. Until then, do not quote absolute probabilities as equal to MITRE’s.

13Limitations

• Flat terrain; no building downwash or wakes.
• One wind direction per hour for the plume path (taken at the mid-height of the trajectory).
• The upset model is simplified pending a fuller Houbolt & Sen treatment: one degree of freedom, strip-theory damping, no pilot input.
• Typical equipment values are placeholders until manufacturer data is supplied.
• NAD83 ≈ WGS84 at the metre level.
• Clusters are treated as independent; separate plumes do not interact.
• Buoyancy uses reduced gravity referenced to the plume temperature, as in Briggs' F = g V r² (T_s − T_a)/T_s.
• Plumes are integrated only while rising. A plume that overshoots and sinks back (a collapsing fountain), and its downdraft, is out of scope.
• Long rows of adjacent units are merged in segments of at most 60 m, each an equivalent round stack, rather than as a true line source.
• Upset gust width is 12.5 c̄ in total, following MITRE's stated assumption. FAR 25.341's legacy gust puts 12.5 c̄ at the peak (25 c̄ total), which would roughly halve the upset threshold. This is pending a check against MITRE output.
• Weather above the top of the profile is held constant (a MITRE CSV stops at 914 m).
• Results depend on the weather model, whose grid cell may differ from the site.

14References

1. D. Gouldey, J. Hopper, J. Schwalbe, September 2012, “Expanded Model for Determining the Effects of Vertical Plumes on Aviation Safety.” MP120461, The MITRE Corporation, McLean, VA.
2. D. Gouldey, J. Hopper, and J. Schwalbe, June 2012. “Validation of the SAIC Model for Predicting the Effects of Vertical Plumes on Aviation Safety.” MP120264, The MITRE Corporation, McLean, VA.
3. K.T. Spillane, “The Rise of Wet Plumes – Conservation Equations and Entrainment Assumptions.” Report No. S0/80/10, State Electricity Commission of Victoria, Research and Development Department, Engineering Research Division, 1980.
4. L.V. Schmidt, 1998. Introduction to Aircraft Flight Dynamics, AIAA Education Series, 397.
5. J.C. Houbolt and A. Sen, 1971. “Single-Degree-of-Freedom Roll Response Due to Two-Dimensional Vertical Gusts.” NASA CR-111966. Aeronautical Research Associates of Princeton, Inc., Princeton, NJ.
6. P. N. Papanicolaou and E. J. List, 1988. “Investigations of round vertical turbulent buoyant jets.” Journal of Fluid Mechanics, vol. 195, 341-391.
7. NOAA, “The Rapid Update Cycle.” http://ruc.noaa.gov/.
8. Cessna, October 2011, “Pilot’s Operating Handbook and Flight Training Supplement Cessna Model 162 – Revision 4.”
9. S. Gudmundsson, 2019. General Aviation Aircraft Design. Butterworth-Heinemann, Waltham, MA.
10. D. Gouldey and J. Hopper, “Determining the Risk of Experiencing Severe Turbulence when Flying through an Exhaust Plume.” International Conference on Research in Air Transportation (ICRAT), Istanbul, Turkey, May 26-30, 2019.
11. G.A. Briggs, 1969. Plume Rise. AEC Critical Review Series, TID-25075, U.S. Atomic Energy Commission.
12. G.A. Briggs, 1984. “Plume rise and buoyancy effects.” In D. Randerson (ed.), Atmospheric Science and Power Production, DOE/TIC-27601, U.S. Department of Energy.
13. D.P. Hoult, J.A. Fay, L.J. Forney, 1969. “A theory of plume rise compared with field observations.” Journal of the Air Pollution Control Association, 19(8), 585–590.
14. B.R. Morton, G.I. Taylor, J.S. Turner, 1956. “Turbulent gravitational convection from maintained and instantaneous sources.” Proceedings of the Royal Society A, 234, 1–23.
15. 14 CFR §23.341, Gust loads factors. U.S. Federal Aviation Regulations.

Back to overview