Back to Odor Forecast

Atmospheric Dispersion Heuristic for Real-Time Odor Estimation: Chemical-Specific Gaussian Plume Implementation with AI-Assisted Development

Authors: Community Environmental Monitoring Project, Alexandria–Pineville, Louisiana
Date: March 2026
Version: 1.4 — Experimental Draft for Peer Review
Experimental Draft for Peer Review — This document describes a heuristic model under active development. Findings are preliminary and should not be cited as established conclusions.

Abstract

This paper presents a novel, web-based atmospheric odor prediction system designed to forecast community-level odor events originating from two wood-treatment industrial facilities in central Louisiana. The system integrates real-time National Weather Service (NWS) forecast data, chemical-specific volatilization physics, Gaussian plume dispersion modeling, terrain-aware transport calculations, and temperature inversion detection into a unified scoring model that produces hourly odor intensity forecasts across a 48-hour horizon. A key innovation is the explicit modeling of open-air treated wood storage yards as dominant area emission sources, with solar-heated surface temperature calculations driving nonlinear volatilization curves that shift source dominance between daytime (area sources) and nighttime (point and elevated sources). The system also incorporates product-specific surface-area-to-volume ratios and diurnally shifting source weights. The application was developed entirely through AI-assisted coding using Claude, Anthropic’s large language model, and has shown high correlation with community-reported odor events as validated through field observations and remote sensing validation of publicly available aerial imagery.

1. Introduction

1.1 Problem Statement

Wood-treatment facilities that use creosote, pentachlorophenol, and other chemical preservatives emit volatile organic compounds (VOCs) that can produce persistent, community-impacting odors. The Alexandria–Pineville metropolitan area in central Louisiana hosts two such facilities that collectively emit approximately 97 tons of VOCs per year. Community members have long reported episodic odor events, but the occurrence and intensity of these events are highly dependent on meteorological conditions, making them difficult to predict.

Traditional air quality monitoring relies on fixed-position sensors and regulatory dispersion models (e.g., AERMOD, CALPUFF) that require professional operation, expensive instrumentation, and are not typically available to the public in real-time. There exists a gap between industrial emissions data publicly available through environmental permits and the community’s ability to anticipate when those emissions will produce noticeable odors at ground level.

1.2 Objective

This project sought to bridge that gap by developing an accessible, browser-based prediction tool that:

  1. Translates publicly available emissions inventory data into real-time odor forecasts
  2. Accounts for the chemical-specific behavior of different VOC compounds
  3. Integrates freely available NWS weather forecast data as the atmospheric driver
  4. Provides spatially resolved predictions on an interactive map
  5. Operates without specialized hardware, software licenses, or atmospheric science training

1.3 AI-Assisted Development

The entire application — encompassing atmospheric physics, chemical engineering principles, geographic information systems, and web development — was developed through iterative collaboration with Claude, Anthropic’s AI assistant. This represents a significant demonstration of how large language models can synthesize domain knowledge across multiple scientific and engineering disciplines to produce functional, accurate software systems. The development process involved translating environmental permit data, atmospheric science principles, and field observations into working algorithms through natural-language dialogue, with the AI generating all source code, mathematical models, and visualization components. The author directed the synthesis of the underlying mathematical models, ensuring that the AI-assisted code implementation correctly reflects established principles of atmospheric physics and chemical volatilization. All physical constants, meteorological logic, and model parameters were reviewed and validated by the project author against published literature to ensure scientific accuracy.


2. Facility Characterization and Emissions Profiles

2.1 Data Sources and Remote Sensing Validation

Emissions data were derived from Louisiana Department of Environmental Quality (LDEQ) permit records, which provide annual emissions estimates by compound for each facility. These permits represent the authoritative public record of expected emissions under normal operating conditions. However, traditional regulatory models often focus on point-source stack emissions, and permit-level fugitive emission estimates may not fully characterize the spatial extent and geometry of ground-level area sources.

To ensure model fidelity, permit data was supplemented by remote sensing validation using publicly available aerial imagery (Google Maps/Earth, 2024–2026) and Street View photography. This analysis enabled the identification and geometric characterization of fugitive area sources — specifically open-air storage yards — which are primary drivers of ground-level odor events but are often underrepresented in stack-height-based dispersion modeling. All imagery analyzed is publicly accessible; source characterizations are derived from observable spatial extent, albedo contrast, and product geometry visible in the imagery. The distinction between area sources (open-air wood storage yards) and point sources (treatment cylinders, kiln stacks) is critical to the model’s emission factor architecture.

Permit-only data policy and conservative bias: To maintain scientific integrity and regulatory traceability, the model intentionally excludes uncharacterized sources. All emission quantities and source types (fugitive vs. stack) used in the scoring algorithm are derived exclusively from LDEQ Electronic Document Management System (EDMS) permit filings. While additional potential area sources (such as unpermitted retention areas or staging surfaces) may be identifiable via remote sensing, they are not included in emission calculations. This ensures that every high-intensity forecast is directly traceable to authorized emissions listed in public regulatory records. Remote sensing validation is used solely for geometric characterization of permitted source categories — not for introducing additional emission sources beyond the permit record.

2.2 Alexandria Facility: Source Characterization and Geometry (LDEQ Permit #2360-00032-07)

The Alexandria facility, located at 72 ft elevation on the Red River valley floor, treats railroad ties with coal-tar creosote using four high-pressure treatment cylinders operating at approximately 140°F. The facility operates continuously (24/7/365) with 7–8 cylinder unloading events per day.

Remote sensing validation of publicly available aerial imagery identifies two distinct emission source categories, characterized for modeling purposes as follows:

Untreated railroad ties are also stored on site in tightly packed, uniform square stacks, primarily in the northern portion of the yard. These lighter-colored stacks are not a significant odor source.

Table 1: Alexandria Facility Emissions Profile

CompoundAnnual Emissions (tons/yr)Fugitive FractionOdor Character
Naphthalene9.2680% (7.41 tons)Mothball/tar
Cresol0.14HighMedicinal/phenolic
Quinoline0.36ModerateChemical
Biphenyl0.20ModeratePleasant/floral
PAHs (mixed)0.26LowTar
Phenol0.09HighSweet/medicinal
Total VOCs18.67
Non-VOC emissions:
PM10 (particulate matter)12.22

The primary odor-surrogate identified for the Alexandria facility is naphthalene, with 80% of permitted emissions categorized as fugitive (originating from equipment interfaces, valves, pumps, and post-treatment product staging rather than from controlled stacks). This ground-level emission profile is a critical variable in the dispersion model, as these releases occur below the atmospheric boundary layer and are immediately influenced by surface-level meteorology. The open-air storage yard is characterized in this model as a large area source with zero effective stack height. This configuration indicates that volatilized naphthalene is subject to horizontal transport even under low-velocity wind conditions, lacking the vertical buoyancy and initial dilution typically provided by elevated stack releases.

2.3 Pineville Facility: Multi-Regime Emission Dynamics (LDEQ Permit #2360-00004-07)

The Pineville facility presents a substantially more complex emissions profile. Located at 103 ft elevation on the Pineville bluff — approximately 30 ft above the Red River valley floor — the facility treats utility poles with multiple preservative systems (creosote, pentachlorophenol, DCOI, and CCA) and operates three treatment cylinders, two dry kilns, and a large tank farm (>103,000 gallons of creosote and >62,000 gallons of diesel carrier solvent).

Remote sensing validation of publicly available aerial imagery indicates a more complex, multi-regime emission profile. Dozens of creosote-treated telephone pole stacks are stored outdoors across several acres of open yard, identifiable in satellite imagery as dark black/brown bundles distinct from the lighter raw/untreated poles also staged on site. These are characterized for modeling purposes as solar-driven area sources; while their cylindrical geometry yields a lower SA:V than railroad ties, their dark albedo (α ≈ 0.05–0.15) drives the solar volatilization model used in Section 4. Kiln stacks and rooftop exhaust vents are also visible, with an estimated effective stack height of 25–32 ft, constituting continuous elevated point sources during operation. The facility’s three source categories — open-air pole yard, ground-level process emissions, and elevated kiln stacks — each exhibit fundamentally different emission and transport behaviors.

Table 2: Pineville Facility Emissions Profile

CompoundAnnual Emissions (tons/yr)Release TypeOdor Character
Kiln Stack Emissions (elevated, 25–32 ft):
Acetaldehyde2.10StackPungent/fruity
Formaldehyde0.86StackSharp/irritating
Methanol0.82StackSweet/pungent
Other kiln VOCs~57.90StackMixed
Ground-Level Emissions:
Cresol1.62FugitiveMedicinal/phenolic
Naphthalene0.84FugitiveMothball
Diesel (carrier)SignificantFugitivePetroleum
Biphenyl0.19FugitivePleasant
Quinoline0.11FugitiveChemical
n-Hexane0.26FugitiveGasoline-like
Total VOCs78.68

The critical insight incorporated into the model is that Pineville’s emissions operate as a three-regime system. Open-air treated pole storage constitutes a solar-driven area source that dominates during daytime hours. Ground-level process emissions (cylinders, tanks, drip pads) contribute cresol, diesel, and naphthalene as fugitive sources. Kiln stacks (25–32 ft) release 78% of total VOCs as light, buoyant gases (molecular weight 30–44 g/mol) that dominate the nighttime and inversion-trapped odor profile. The treated poles have a lower surface-area-to-volume ratio than Alexandria’s railroad ties (cylindrical vs. flat rectangular geometry), resulting in a lower per-acre area source emission rate despite covering a comparable footprint. These three regimes exhibit fundamentally different emission rates, temperature sensitivities, and atmospheric transport behavior, requiring separate modeling approaches with diurnally shifting source weights.


3. Meteorological Data Integration

3.1 Data Sources and Grid Configuration

The system integrates three complementary weather data sources: NWS hourly forecast data, Open-Meteo model forecasts for multi-model blending and gap-filling, and real-time KAEX observation data for current-hour accuracy. The NWS provides hourly forecast data on a 2.5 km grid covering the continental United States. The system queries the grid point nearest to the Alexandria–Pineville area:

3.2 Forecast Parameters

The following parameters are extracted from NWS hourly forecasts for each prediction period:

ParameterSource EndpointResolution
Temperature (°F)/gridpoints/LCH/88,143/forecast/hourlyHourly
Dewpoint (°F)/gridpoints/LCH/88,143 (gridData)Hourly
Relative Humidity (%)/gridpoints/LCH/88,143/forecast/hourlyHourly
Wind Speed (mph)/gridpoints/LCH/88,143/forecast/hourlyHourly
Wind Direction (compass)/gridpoints/LCH/88,143/forecast/hourlyHourly
Short Forecast (text)/gridpoints/LCH/88,143/forecast/hourlyHourly
Area Forecast Discussion/products/types/AFD/locations/LCH~6 hourly

Wind speed values, which are reported in range format (e.g., “5 to 10 mph”), are parsed by extracting all numeric values and using the maximum, providing a conservative estimate for dilution calculations.

3.2.1 Observation Blending

For the current hour, the system replaces NWS forecast values with actual measurements from the KAEX ASOS (Automated Surface Observing System) station. This “observation blending” addresses a fundamental limitation of forecast-only models: NWS forecasts are predictions that may diverge significantly from actual conditions, particularly for wind speed and direction. In initial testing, forecast wind speed errors of 5+ mph were observed, which can substantially alter plume dispersion calculations.

The following additional parameters are available from KAEX observations but not from NWS forecasts:

ParameterSourceUse in Model
Barometric Pressure (mb)/stations/KAEX/observations/latestInversion detection: pressure trend signal
Visibility (miles)/stations/KAEX/observations/latestInversion detection: direct trapping measurement

For all past and current hours in the 48-hour timeline, the system fetches the full range of archived KAEX observations from the D1 database and blends each period with its closest observation (within a 1-hour tolerance window). This means scrolling back through the timeline shows actual measured conditions rather than stale forecast predictions. Each blended period is flagged (hasObservation = true) to enable observation-aware weighting in downstream algorithms (see Section 4.4) and to indicate the data source in the UI (e.g., “Current Conditions (Observed)” vs. “Forecast at 9 PM”).

3.2.2 Near-Term Forecast Bias Correction

NWS forecasts can diverge significantly from actual conditions, particularly during transitional weather. In field testing on March 17, 2026, the overnight NWS forecast predicted calm winds and conditions favorable for an odor event. The model generated an elevated odor score, but by morning the observed winds at KAEX were substantially different — the predicted event did not materialize. The forecast had been wrong for hours, yet the model continued to predict elevated odor for the next several hours based on the same stale forecast.

To address this, the system implements observation-based bias correction for near-future forecast periods. When the current hour has both an observation and the original forecast values, the system computes the forecast error (delta) for each parameter:

$$\Delta_p = \text{Observed}_p - \text{Forecast}_p$$

where $p \in \{\text{temperature, wind speed, wind direction, humidity}\}$. This delta is then applied to the next 2 forecast hours with exponentially decaying weight:

Hours AheadCorrection WeightRationale
+1 hour75%Strong persistence — conditions rarely reverse within an hour
+2 hours33%Moderate persistence — errors begin to self-correct
+3 hours and beyond0%Raw forecast — bias may no longer apply

For scalar parameters (temperature, wind speed, humidity), correction is arithmetic:

$$\text{Corrected}_p = \text{Forecast}_p + \Delta_p \times w$$

Wind direction requires angular interpolation to handle compass wraparound (e.g., a delta between 350° and 10° is 20°, not −340°). The delta is normalized to the [−180°, +180°] range before applying the weighted correction, and the compass direction label is recalculated from the corrected degrees.

Wind speed is clamped to a minimum of 0 mph, and humidity is clamped to [0%, 100%] after correction. Bias-corrected periods are flagged (hasBiasCorrection = true) and displayed in the UI as “Forecast at 9 PM (Corrected)” to distinguish them from raw NWS forecasts without cluttering the interface. The correction deltas are logged to the browser console for diagnostic purposes.

This approach is conservative by design. Two hours of correction with rapid decay means the system quickly adapts when reality diverges from the forecast, but does not over-correct for distant hours where the NWS model may have self-corrected. The result is that the “shoulder hours” immediately ahead of the current time — the hours most relevant to a user deciding whether to close their windows — reflect the actual atmospheric trajectory rather than a potentially stale prediction.

3.2.3 Multi-Model Forecast Blending (NWS + Open-Meteo)

NWS hourly forecasts occasionally exhibit gaps or stale data, particularly during the critical 6–8 AM window when surface inversions are forming or breaking and small forecast errors can produce large prediction errors in plume behavior. To improve forecast coverage and accuracy, the system supplements NWS data with forecasts from the Open-Meteo API, a free service that aggregates output from global numerical weather prediction models including ECMWF IFS and NCEP GFS.

Both forecast sources are fetched in parallel during each 30-minute cron cycle and stored separately in D1, preserving raw data from each source for accuracy comparison over time. Blending occurs at serve time when the frontend requests forecast data, allowing blend weights to be tuned without re-fetching historical data. The blending logic time-aligns periods from both sources by UTC hour key and applies the following rules:

ConditionResultSource Tag
NWS period exists, Open-Meteo missingUse NWS as-is"nws"
Open-Meteo exists, NWS missingUse Open-Meteo (fills gap)"openmeteo"
Both sources availableTime-dependent weighted average (see below)"blended"

Time-Dependent Blend Weights

Rather than applying a fixed blend ratio, the system uses time-dependent weighting that shifts from NWS-heavy at short lead times to Open-Meteo/ECMWF-heavy at longer lead times. This exploits the complementary strengths of each source: NWS forecasts incorporate local KAEX station data and regional forecaster adjustments from the Lake Charles (LCH) and Shreveport (SHV) offices, providing superior mesoscale accuracy for the first 6–12 hours. Beyond 24 hours, however, ECMWF IFS—the primary model behind Open-Meteo—consistently outperforms GFS-derived NWS forecasts in synoptic-scale verification scores for temperature and wind.

The NWS weight wNWS is computed as a function of forecast lead hours h:

$$w_{\text{NWS}}(h) = \begin{cases} 0.70 & h \leq 6 \\ 0.70 - \dfrac{0.40}{42}(h - 6) & 6 < h < 48 \\ 0.30 & h \geq 48 \end{cases}$$

This produces a linear ramp from 70% NWS at 0–6 hours to 30% NWS at 48 hours:

Lead TimeNWS WeightOpen-Meteo Weight
0–6 h70%30%
12 h64%36%
24 h53%47%
36 h41%59%
48 h+30%70%

The 6-hour flat zone at the start avoids unnecessary churn when both sources agree closely in the very near term. Each blended period includes a nwsWeight field in the API response, making the per-hour weighting transparent to the frontend and useful for debugging.

Blending Formulas

For numeric parameters (temperature, dewpoint, humidity, wind speed), blending is arithmetic using the time-dependent weight:

$$\text{Blended}_p = w_{\text{NWS}} \times \text{NWS}_p + (1 - w_{\text{NWS}}) \times \text{OM}_p$$

Wind direction blending uses vector averaging to handle compass wraparound correctly:

$$\theta_{\text{blend}} = \text{atan2}\!\left(w_{\text{NWS}} \sin\theta_{\text{NWS}} + w_{\text{OM}} \sin\theta_{\text{OM}},\; w_{\text{NWS}} \cos\theta_{\text{NWS}} + w_{\text{OM}} \cos\theta_{\text{OM}}\right)$$

Non-numeric fields (shortForecast, icon) always come from NWS, as Open-Meteo does not provide text descriptions. Open-Meteo contributes an additional cloudCover field (0–100%) to forecast periods, which is more granular than the NWS text-based sky condition and may be used to improve future-hour inversion detection.

The Open-Meteo API is requested with timeformat=unixtime to receive Unix timestamps, avoiding timezone ambiguity in the Cloudflare Worker’s UTC runtime environment. Timestamps are converted to ISO 8601 format for consistent alignment with NWS periods, which include timezone offsets. The API response includes a sources metadata object with nws_fetched_at and openmeteo_fetched_at timestamps, enabling the frontend to verify data freshness from both sources.

3.3 D1-Backed Data Architecture

Early versions of the system used a Cloudflare Worker with KV (key-value) storage as a caching proxy for the NWS API. Each user page load could trigger NWS API requests through this proxy, creating rate-limit risk and coupling frontend availability to NWS uptime. This architecture was replaced with a fully D1-backed design where no user action ever triggers an NWS API call.

The current architecture uses a two-tier strategy:

  1. Browser-side cache: Hourly forecasts are cached for 30 minutes, Area Forecast Discussions for 1 hour, and observation data for 10 minutes, using in-memory storage. This reduces repeated API calls within a single browser session.
  2. D1 database (Cloudflare Workers + Cron): A scheduled cron trigger runs every 30 minutes and is the sole consumer of external weather APIs. Each cycle fetches four data sources in parallel and archives them to D1:
    • KAEX observations: Latest surface measurements (temperature, dewpoint, humidity, wind, pressure, visibility) stored alongside the corresponding NWS forecast period for the same time, enabling forecast accuracy tracking.
    • NWS hourly forecast: The full 48–156 period NWS hourly forecast is pre-parsed (unit conversions, wind speed parsing, compass-to-degrees lookup) and stored as a JSON blob in the hourly_forecasts table.
    • Open-Meteo forecast: A 48-hour forecast from the Open-Meteo API (ECMWF/GFS models) is parsed from parallel arrays into the same period shape as NWS data and stored in a separate openmeteo_forecasts table. Both sources are preserved independently for accuracy comparison; blending occurs at serve time (see Section 3.2.3).
    • Area Forecast Discussion: The latest AFD text is archived with deduplication by NWS product ID, so repeated cron cycles don’t re-fetch unchanged discussions.

The frontend reads exclusively from /api/* endpoints backed by D1. Response times are sub-millisecond for D1 reads, compared to 200–2000 ms for live NWS API calls. The system is resilient to NWS outages: if the cron cycle fails, the previous forecast remains available in D1 until the next successful fetch.

Automatic cleanup keeps storage bounded: NWS and Open-Meteo forecast snapshots older than 7 days and AFD products older than 30 days are pruned each cron cycle. Observation data is retained indefinitely for long-term analysis. At approximately 200 bytes per observation and 48 observations per day, annual storage is approximately 1.7 MB — negligible cost for a dataset that becomes increasingly valuable for pattern analysis, seasonal calibration, and forecast error characterization.


4. Core Prediction Model

4.1 Multiplicative Factor Architecture

The system computes odor intensity as a dimensionless score (0–100) using a multiplicative combination of independent environmental factors:

$$\text{Score} = E \times T \times I \times D \times H \times G \times 100$$

Where:

Each factor is normalized to a baseline of approximately 1.0 under typical conditions, with values below 1.0 indicating conditions that reduce odor impact and values above 1.0 indicating amplification. The multiplicative structure ensures that any single factor near zero (e.g., strong wind providing high dilution) can suppress the overall score regardless of other favorable conditions — reflecting the physical reality that effective dispersion through any mechanism reduces ground-level concentrations.

4.2 Emission Factor (E)

The emission factor models the rate at which VOCs volatilize from facility surfaces, equipment, and liquid impoundments as a function of temperature and time of day. The model distinguishes between two source categories at Alexandria — area sources (open-air treated wood storage) and point sources (treatment cylinders) — and three at Pineville, which adds elevated sources (kiln stacks). Source weighting shifts diurnally: solar-heated area sources dominate during daytime hours, while point and elevated sources dominate at night and during inversions.

4.2.1 Solar Surface Temperature Model

Dark, creosote-saturated surfaces absorb 85–95% of incident solar radiation (albedo 0.05–0.15), consistent with published values for dark organic coatings and asphalt. The resulting surface temperature exceeds ambient air temperature by a margin that increases nonlinearly with ambient temperature, consistent with the Stefan-Boltzmann law (radiative heating scales as $T^4$, while convective cooling scales linearly).

Solar Position Calculation

Rather than using fixed daytime hours (e.g., 8 AM–6 PM), the model calculates astronomically correct sunrise and sunset times for each forecast period using the solar declination equation (Spencer, 1971). This accounts for seasonal variation in day length and DST transitions:

The solar declination $\delta$ is computed from the day of year $n$:

$$\gamma = \frac{2\pi(n-1)}{365}$$ $$\delta = 0.006918 - 0.399912\cos\gamma + 0.070257\sin\gamma - 0.006758\cos 2\gamma + \ldots$$

The hour angle at sunrise/sunset is:

$$\omega_0 = \arccos(-\tan\phi \cdot \tan\delta)$$

where $\phi$ is the site latitude (31.31°N). Solar noon is corrected using the Equation of Time, and sunrise/sunset times are derived from solar noon ± $\omega_0 / 15$ hours.

Solar Elevation Scaling

The solar heating boost is not applied as a flat on/off switch but is scaled by the normalized solar elevation angle $\alpha$, which represents the sun’s angular height above the horizon relative to its peak for the day:

$$S_I = \frac{\alpha(h)}{\alpha_{\max}}$$

where $\alpha(h)$ is the solar elevation at hour $h$ and $\alpha_{\max}$ is the peak elevation at solar noon. This ensures the heating boost ramps gradually after sunrise (low sun angle = weak heating), peaks at solar noon, and tapers toward sunset. The effective surface temperature becomes:

$$T_{\text{eff}} = T_F + \Delta T_{\text{solar}} \cdot S_I$$

where $\Delta T_{\text{solar}}$ is the peak solar heating boost (at full intensity), which varies by product type:

Table 3: Peak Solar Heating Boost by Ambient Temperature and Product Type

Ambient TempRailroad Ties (Alexandria)Telephone Poles (Pineville)Rationale
> 85°F40–62°F boost30–48°F boostSharp nonlinear ramp; surfaces reach 125–155°F at solar noon
70–85°F15–40°F boost10–30°F boostModerate solar heating
55–70°F0–15°F boost0–10°F boostMild solar effect
< 55°F0°F0°FInsufficient solar input

These peak values are attained at solar noon ($S_I = 1.0$). At one hour after sunrise, $S_I$ is typically 0.2–0.3, so a 95°F day would produce only ~12°F of surface boost rather than the full 55°F. This prevents the model from applying unrealistic solar heating during early morning or late evening when the sun is near the horizon.

Railroad ties receive a larger solar boost than telephone poles due to their higher surface-area-to-volume ratio. A tie (approximately 7” × 9” × 8.5’) exposes roughly 2–3× more creosote-saturated surface per unit volume than a cylindrical pole (approximately 12” diameter × 40’). More exposed surface area means more evaporative surface, and the flat horizontal faces of stacked ties receive more direct-normal solar radiation than the curved surfaces of cylindrical poles. At night ($S_I = 0$), surfaces cool to near ambient temperature and $T_{\text{eff}} = T_F$.

4.2.2 Alexandria (Two-Source Model)

The Alexandria emission model combines two source categories with diurnally shifting weights:

Source 1: Open-air tie storage (area source). Vapor pressure is calculated using the Clausius-Clapeyron relation at the effective surface temperature:

$$f_{\text{ties}} = \text{clamp}\!\left(2^{(T_{\text{eff}} - 77) / 15},\; 0.35,\; 5.0\right)$$

The floor of 0.35 reflects field observations that creosote-soaked wood emits detectable naphthalene even at 40°F. The ceiling is raised to 5.0 (from the previous 3.0) to accommodate the high effective surface temperatures driven by solar heating on hot summer days.

Source 2: Cylinder operations (point source). Internal temperature of 140°F is maintained regardless of ambient conditions. The factor is constant at 1.0, reflecting the 24/7 nature of the 7–8 daily unloading events.

The combined Alexandria emission factor uses diurnally shifting source weights:

$$E_{\text{Alex}} = B_I \times (f_{\text{ties}} \times w_1 + f_{\text{cyl}} \times w_2)$$

Table 4: Alexandria Source Weights by Time of Day

SourceDaytime (sunrise–sunset)NighttimeRationale
Open-air tie yard0.650.40Solar heating drives daytime off-gassing; cools at night but residual emission persists
Cylinder operations0.350.60Constant 24/7; relatively more important when area sources cool

where $B_I = 0.5$. Both facilities operate 24/7/365 with no publicly known shift or batch schedule, so the model does not assume peak operational hours. Temporal variation in odor intensity is instead driven entirely by atmospheric factors: temperature-dependent volatilization, inversion trapping, solar heating, and the diurnal factor (Section 4.5).

4.2.3 Pineville (Three-Source Model)

The Pineville emission model separately calculates three source contributions with diurnally shifting weights:

Source 1: Open-air pole storage (area source). Vapor pressure is calculated at effective surface temperature using the cresol curve (doubles every 18°F) since cresol dominates the odor profile of treated poles due to its extraordinarily low detection threshold (~0.001 ppm):

$$f_{\text{poles}} = \text{clamp}\!\left(2^{(T_{\text{eff}} - 77) / 18},\; 0.25,\; 4.0\right)$$

Source 2: Ground-level process emissions (point/fugitive sources). This regime models emissions from cylinders, storage tanks, drip pads, and fugitive sources using a weighted blend of compound-specific factors:

A tank breathing modifier adds 20% during afternoon hours (1–5 PM) when temperatures exceed 70°F, modeling thermal expansion of liquid in the facility’s 165,000+ gallon tank farm. Cylinder operations contribute a constant baseline blended at 25% weight.

Source 3: Kiln stack emissions (elevated source, 25–32 ft). Uses a base factor of 0.7 reflecting reduced ground-level impact of elevated releases, with adjustments:

The combined Pineville emission factor uses diurnally shifting source weights:

$$E_{\text{Pine}} = B_I \times (f_{\text{poles}} \times w_1 + f_{\text{ground}} \times w_2 + f_{\text{kiln}} \times w_3)$$

Table 5: Pineville Source Weights by Time of Day

SourceDaytime (sunrise–sunset)NighttimeRationale
Open-air pole yard0.300.10Major daytime source via solar heating; cools at night
Ground-level process0.250.20Continuous but temperature-modulated
Kiln stacks0.450.70Dominates by mass (78% of VOCs); elevated emissions get trapped during inversions

where $B_I = 0.55$. The kiln retains the highest weight even during daytime because it accounts for 78% of total facility VOC emissions by mass, but the open-air pole yard’s weight increases during daytime to reflect the solar-driven amplification of area source off-gassing.

4.3 Transport Factor (T)

The transport factor models how emissions move through the atmosphere from the source to a receptor point, using a simplified Gaussian plume approach.

4.3.1 Wind-Driven Transport (Wind Speed > 3 mph)

When definable wind exists, the transport factor combines three components:

Angular dispersion uses a Gaussian distribution centered on the downwind direction:

$$f_{\text{angular}} = \exp\!\left(-\frac{1}{2}\left(\frac{\Delta\theta}{\sigma}\right)^{\!2}\right)$$

where $\Delta\theta$ is the angular offset between the plant-to-receptor bearing and the downwind direction, and $\sigma$ is the plume spread parameter. The spread parameter varies with wind speed to reflect turbulence effects:

Wind Speedσ (degrees)Max Transport Distance
≤ 8 mph50°5 km
8–15 mph30°7 km
> 15 mph20°10 km

At low wind speeds, greater atmospheric meandering produces wider plume spread (larger σ), while stronger winds create more directionally focused transport with greater reach.

Distance decay follows a power-law relationship whose exponent varies with atmospheric stability (see Section 4.3.3):

$$f_{\text{distance}} = \left(1 - \frac{d}{d_{\max} \times M_s}\right)^{\!\beta(s)}$$

where $d$ is the Haversine distance from the plant to the receptor point, $d_{\max}$ is the wind-speed-dependent maximum transport distance, $M_s$ is a stability-dependent range multiplier, and $\beta(s)$ is the stability-dependent decay exponent. Under neutral conditions ($s = 0.5$), the exponent reduces to the original field-calibrated value of 1.8.

Wind dilution models the turbulent mixing effect of wind, modulated by atmospheric stability:

$$f_{\text{dilution}} = 1 - \left(1 - \frac{5}{v_w}\right) \times \eta(s)$$

where $v_w$ is wind speed in mph and $\eta(s) = 1.0 - 0.35s$ is the stability-dependent dilution efficiency. In unstable conditions ($\eta \approx 1$), wind produces full turbulent dilution. In stable conditions ($\eta \approx 0.65$), vertical mixing is suppressed by the BVF-limited $\sigma_z$ ceiling (Section 4.3.3), so wind primarily advects the plume along its trajectory rather than dispersing it vertically.

The complete wind-driven transport factor is:

$$T = f_{\text{angular}} \times f_{\text{distance}} \times f_{\text{dilution}}$$

4.3.2 Calm-Air Transport (Wind Speed ≤ 3 mph)

Under calm conditions, the directional plume model is replaced with an isotropic (circular) pooling pattern. Emissions spread uniformly in all directions with a stability-adjusted maximum radius and decay exponent:

$$f_{\text{distance, calm}} = \left(1 - \frac{d}{4.0 \times M_s}\right)^{\!\max(0.8,\; \beta(s) - 0.4)}$$

The additional exponent reduction of 0.4 (floored at 0.8) reflects the compounding effects of calm conditions: no wind-driven mixing and no directional transport, so the pooling pattern decays even more gently than a wind-driven plume at the same stability level. During strong inversions ($s \approx 0.95$), the calm radius extends to approximately 5.6 km with a decay exponent near 0.8, consistent with the 2026-03-13 field observation of odor detection at ~3 km under calm/inversion conditions.

4.3.3 Stability-Dependent Dispersion (Brunt-Väisälä Vertical Confinement)

In the EPA’s AERMOD regulatory dispersion model, the vertical dispersion parameter $\sigma_z$ does not grow without bound. In stable atmospheric conditions, the Brunt-Väisälä frequency (BVF) — a measure of the atmosphere’s resistance to vertical displacement — imposes an asymptotic ceiling on vertical plume spread:

$$N = \sqrt{\frac{g}{T}\frac{\partial\theta}{\partial z}}$$ $$\sigma_z \;\xrightarrow{x \to \infty}\; \sqrt{\frac{0.54\,\sigma_w}{N}}$$

where $N$ is the BVF, $g$ is gravitational acceleration, $T$ is absolute temperature, $\partial\theta/\partial z$ is the potential temperature gradient (positive in stable conditions), and $\sigma_w$ is the vertical velocity variance. When $N$ is large (strong stability), $\sigma_z$ asymptotes to a small value, confining the plume to a thin vertical layer. When the atmosphere is unstable ($\partial\theta/\partial z < 0$), $N$ is undefined (imaginary), and $\sigma_z$ grows freely with downwind distance.

The physical consequence is that ground-level concentration $C$ at large distances behaves differently under stable vs. unstable conditions. In the standard Gaussian plume model:

$$C \propto \frac{1}{\sigma_y \cdot \sigma_z}$$

Under unstable conditions, both $\sigma_y$ and $\sigma_z$ grow with distance, producing rapid $\sim 1/x^2$ concentration decay. Under stable conditions, $\sigma_z$ is capped by BVF while $\sigma_y$ continues to grow, producing slower $\sim 1/x$ decay — the plume stays concentrated and detectable at significantly greater distances. This is the fundamental physics behind why temperature inversions cause odor complaints far from the source.

This system does not have the meteorological profile data (vertical temperature soundings, friction velocity $u_*$, Monin-Obukhov length $L$) required to compute BVF directly. Instead, we estimate an atmospheric stability parameter $s \in [0, 1]$ from available NWS signals and use it to parameterize the BVF effect on dispersion:

Stability estimation. The stability parameter integrates three signal categories:

  1. Inversion score (primary): The multi-signal inversion detection system (Section 4.4) provides a composite stability proxy. A high inversion score implies strong stratification (high BVF): $$s_{\text{base}} = 0.5 + 0.45 \times S_{\text{inv}}$$
  2. Solar heating (daytime cap): Convective instability during hot afternoons limits stability regardless of other signals. When $T_F > 85$°F during peak daylight hours, $s$ is capped at 0.25; when $T_F > 70$°F, $s$ is capped at 0.4.
  3. Mechanical turbulence (wind override): Strong wind generates turbulent kinetic energy that overwhelms thermal stratification. At $v_w \geq 15$ mph, $s$ is capped at 0.35; at $v_w \geq 10$ mph, at 0.45. At night, $s$ has a floor of 0.55 reflecting minimum radiative-cooling stability.

Decay exponent. The stability parameter modulates the distance decay exponent continuously:

$$\beta(s) = 2.2 - 1.2\,s$$

This maps stability to physically motivated decay regimes:

Stability$s$$\beta$Physical Regime
Strongly unstable0.02.2Deep convective mixing; $\sigma_z$ grows rapidly; fast 3D dilution
Neutral0.51.8Moderate mixing; original field-calibrated baseline
Moderately stable0.71.36BVF limits vertical spread; plume concentrated at ground level
Strongly stable0.951.06$\sigma_z$ nearly constant; $\sim 1/x$ decay; maximum range

Range multiplier. The maximum transport distance is scaled by stability to reflect the extended detectability of vertically confined plumes:

$$M_s = 0.85 + 0.55\,s$$
Stability$M_s$Effect on $d_{\max}$
Strongly unstable ($s = 0$)0.8515% shorter range (rapid vertical dilution)
Neutral ($s = 0.5$)1.12~Original range
Strongly stable ($s = 0.95$)1.3737% extended range (BVF vertical confinement)

Together, the reduced decay exponent and extended range under stable conditions produce a compounding effect: concentrations are both higher at any given distance and detectable at greater distances. For example, at 3 km downwind with $d_{\max} = 5$ km and light wind:

This 18× ratio between unstable and stable concentrations at 3 km is consistent with the well-documented observation that industrial odor complaints increase dramatically during nighttime inversions and are rarely reported during hot, windy afternoons — a pattern that the previous fixed-exponent model could only partially capture through the separate inversion multiplier.

4.3.4 Briggs Plume Rise for Kiln Stack

The Pineville facility operates two dry kilns that exhaust through stacks at 25–32 ft (approximately 9 m). Hot exhaust gases are buoyant and rise above the physical stack height before reaching equilibrium with the surrounding atmosphere. The effective release height $H_e = H_s + \Delta H$ determines how much of the plume reaches ground level. This system implements a simplified version of the Briggs plume rise equations used in EPA's AERMOD model (subroutine prise.f).

Buoyancy flux. The thermal buoyancy of the stack exhaust is quantified by:

$$F_b = g \cdot V_s \cdot \left(\frac{D_s}{2}\right)^2 \cdot \frac{T_s - T_a}{T_s}$$

where $g = 9.81$ m/s², $V_s$ is the stack exit velocity (7.5 m/s), $D_s$ is the stack diameter (1.0 m), $T_s$ is the stack exit temperature (~473 K for 200°C kiln exhaust), and $T_a$ is the ambient temperature in Kelvin.

Stable rise (Briggs, 1984):

$$\Delta H_{\text{stable}} = 2.66 \left(\frac{F_b}{U \cdot N^2}\right)^{1/3}$$

where $U$ is the wind speed at stack height (m/s) and $N$ is the Brunt-Väisälä frequency estimated from the stability parameter: $N = \max(0.005,\; 0.04 \times s)$ s¹.

Unstable/neutral rise (Briggs):

$$\Delta H_{\text{unstable}} = \frac{1.6 \cdot F_b^{1/3} \cdot x_f^{2/3}}{U}$$

where $x_f = 49 \cdot F_b^{5/8}$ for $F_b < 55$ m&sup4;/s³. The final rise is blended between stable and unstable using the stability parameter: $\Delta H = s \cdot \Delta H_s + (1-s) \cdot \Delta H_u$, capped at 50 m.

Key behavior:

The Briggs rise replaces the previous fixed ground-level fraction (0.35) for elevated stack releases with a physically-based Gaussian ground-contact formula:

$$f_{\text{ground}} = \exp\!\left(-\frac{H_e^2}{2\sigma_z^2}\right)$$

where $\sigma_z$ is estimated at a typical community distance (~500 m) using stability-dependent growth: $\sigma_z = 8 + 17(1-s)$ meters, ranging from 8 m (stable, BVF-capped) to 25 m (unstable, deep mixing).

4.3.5 Log-Law Wind Profile Adjustment

The KAEX ASOS station measures wind speed at 10 m above ground. The plume at effective height $H_e$ (typically 15–35 m after Briggs rise) experiences a different wind speed due to the atmospheric boundary layer velocity profile. Following AERMOD's wind profile module (windgrid.f), this system applies the neutral Monin-Obukhov log-law:

$$U(z) = U_{\text{ref}} \cdot \frac{\ln(z/z_0)}{\ln(z_{\text{ref}}/z_0)}$$

where $z_0 = 0.5$ m is the surface roughness length (suburban/light industrial land use) and $z_{\text{ref}} = 10$ m is the ASOS measurement height. Representative wind speed ratios:

Height$U(z)/U_{10}$Application
1.7 m (breathing)0.41Ground-level exposure
9 m (stack height)0.96Stack exit conditions
20 m (moderate rise)1.23Plume after modest Briggs rise
30 m (large rise)1.37Plume after significant Briggs rise

The log-law adjustment is applied in the wind dilution calculation: the plume at elevation experiences ~15–35% stronger wind than measured at 10 m, producing additional dilution that reduces ground-level concentrations. The ratio is capped at 2.0 to prevent extreme extrapolation at heights well above the surface layer.

4.3.6 Continuous Meander Transition

At low wind speeds, turbulent eddies are comparable in magnitude to the mean wind, causing the plume to meander randomly in all directions rather than following a coherent downwind path. AERMOD models this transition continuously via the meander fraction (subroutine MEANDR in calc2.f):

$$F_{\text{meander}} = 1 - \text{FRAN}$$

where FRAN (the coherent fraction) is:

$$\text{FRAN} = \frac{2\sigma_v^2 + U_{\text{mean}}^2 \left(1 - e^{-t_{\text{trav}}/T}\right)}{U^2}$$

with $\sigma_v \approx 0.5 \cdot \min(U, 2)$ m/s (lateral turbulence), $U_{\text{mean}}^2 = \max(0.01, U^2 - 2\sigma_v^2)$, $t_{\text{trav}} = d/U$ (travel time), and $T = 86400$ s (24-hour meander timescale, AERMOD default).

The previous model used a hard cutoff at 3 mph: below this threshold, emissions were modeled as pooling circularly; above it, as a directional Gaussian plume. This produced an abrupt behavioral change in the 3–8 mph transition zone. The continuous meander replaces this with a smooth blend:

$$S_{\text{transport}} = F_m \cdot S_{\text{circular}} + (1 - F_m) \cdot S_{\text{directional}}$$

Representative meander fractions:

Wind SpeedDistance$F_m$Behavior
2 mph2 km~0.85Mostly circular pooling
5 mph2 km~0.45Even blend of circular and directional
8 mph3 km~0.15Mostly directional, slight spread
12+ mphany~0.02Nearly pure directional transport

The plume visualization geometry similarly blends between a circular radius and a directional fan: the Gaussian angular spread parameter $\sigma_\theta$ transitions smoothly from its wind-regime value toward 180° (uniform circle) as meander increases: $\sigma_{\theta,\text{eff}} = \sigma_\theta + F_m \cdot (180 - \sigma_\theta)$.

4.3.7 Haversine Distance and Bearing Calculations

All geographic calculations use the Haversine formula for great-circle distance on a spherical Earth:

$$d = 2R \arcsin\!\left(\sqrt{\sin^2\!\left(\frac{\Delta\phi}{2}\right) + \cos\phi_1 \cos\phi_2 \sin^2\!\left(\frac{\Delta\lambda}{2}\right)}\right)$$

where $R = 6{,}371$ km (Earth’s mean radius), $\phi$ denotes latitude, and $\lambda$ denotes longitude. Bearing calculations use the standard forward azimuth formula to determine the direction from plant to receptor.

4.4 Temperature Inversion Detection and Factor (I)

Temperature inversions — stable atmospheric layers where temperature increases with altitude — are the single most important amplifier of ground-level odor concentrations. Under inversion conditions, vertical mixing is suppressed, trapping emissions near the surface and allowing concentrations to build over time.

4.4.1 Multi-Signal Detection

Direct measurement of vertical temperature profiles requires radiosonde (weather balloon) data, which is available only twice daily from stations often hundreds of kilometers away. The system instead detects inversions through nine meteorological proxy signals, each providing independent evidence of atmospheric stability. The signals are designed to capture two distinct inversion mechanisms: fog-type inversions (detected primarily by dewpoint convergence and visibility) and radiation inversions (detected by sky condition and persistence signals), which form under opposite moisture conditions but produce similar ground-level trapping:

Signal 1: Dewpoint-Temperature Convergence (Weight: 0.10)

The spread between temperature and dewpoint indicates moisture saturation proximity. In Louisiana’s humid climate, a baseline nighttime spread of 2–5°F is typical; convergence below this range strongly suggests a stable surface layer with fog-type inversion conditions. Note that this signal’s weight was reduced from its original 0.25 after field observations revealed that radiation inversions — which form under clear, dry conditions with large T−Td spreads — were being missed entirely (see Signal 8):

SpreadSignal ValueInterpretation
≤ 1°F1.0Fog-like, strong inversion
1–2°F0.7Very close, likely inversion
2–4°F0.3Normal night, not indicative
4–8°F0.1Moderately dry
> 8°F0.0Dry, inversions unlikely

Signal 2: Temperature Trend Analysis (Weight: 0.10)

The system tracks temperature changes across forecast periods to identify trend breaks characteristic of inversion formation. During nighttime hours (midnight–7 AM), rising or flat temperatures indicate an inversion layer is insulating the surface from further radiative cooling (signal value 0.9). During post-sunset hours (6–11 PM), a rapid cooling trend followed by a plateau suggests inversion formation (signal value 0.7). A minimum of three prior forecast periods is required for trend detection reliability.

Signal 3: Wind Speed (Weight: 0.20 forecast / 0.15 observation, with VETO capability)

Wind is the strongest single predictor of inversion persistence, as mechanical turbulence from wind disrupts stable atmospheric layers:

Wind SpeedSignal Value
≤ 2 mph1.0
2–5 mph0.7
5–8 mph0.3
≥ 8 mphVETO: total score capped at 0.2

The veto mechanism prevents the system from predicting inversions during windy conditions regardless of other signals, avoiding false positives that might occur on windy, humid nights where dewpoint convergence alone might suggest stability.

Signal 4: Time-of-Day Prior (Weight: 0.05)

A Bayesian prior reflecting the climatological probability of inversions by time of day:

Time PeriodSignal ValueBasis
3–8 AM0.9Maximum radiative cooling, strongest inversions
Midnight–2 AM0.7Stable nocturnal boundary layer
7–11 PM0.6Cooling transition, inversions forming
11 AM–3 PM0.05Peak solar heating, inversions broken

Signal 5: NWS Area Forecast Discussion Keywords (Weight: 0.15 forecast / 0.10 observation)

The system performs natural language parsing of the NWS Area Forecast Discussion (AFD), a free-text narrative written by meteorologists that often describes atmospheric stability conditions not captured in gridded forecast data. The following keywords are scanned:

inversion, fog, stagnation, mixing height, boundary layer, stable, decoupled, trapped, poor dispersion, low-level, nocturnal

Each keyword match contributes 0.25 to the signal value (maximum 1.0), providing a mechanism to incorporate expert meteorological judgment into the automated prediction.

Signal 6: KAEX Visibility (Observation-only, Weight: 0.15)

When KAEX observation data is available, the ASOS visibility sensor provides a direct measurement of atmospheric trapping. Low visibility indicates that particles and droplets cannot disperse vertically — a physical confirmation of inversion conditions that the dewpoint signal can only infer:

VisibilitySignal ValueInterpretation
≤ 0.5 mi1.0Dense fog — strong inversion confirmed
0.5–1 mi0.9Fog
1–3 mi0.7Mist/haze — significant trapping
3–5 mi0.4Light haze — some trapping
5–7 mi0.15Slight reduction — marginal
> 7 mi0.0Clear, good vertical mixing

Signal 7: Barometric Pressure Trend (Observation-only, Weight: 0.10)

Steady or rising barometric pressure with calm winds indicates a stagnant air mass — the synoptic-scale pattern most conducive to persistent inversions. Falling pressure typically indicates an approaching front that will introduce mechanical mixing and break any existing inversion. The signal requires pressure data from at least two observation periods to compute a trend rate:

Pressure Rate (mb/hr)Signal ValueInterpretation
≥ 0.5 (rising fast)0.3High pressure building, moderate stagnation signal
−0.1 to 0.5 (steady)0.7Classic stagnation — air mass not moving
−0.5 to −0.1 (slowly falling)0.3Marginal — change approaching
< −0.5 (falling fast)0.1Front approaching — mixing likely

Signal 8: Sky Condition — Radiation Inversion Detection (Weight: 0.20 forecast / 0.10 observation)

Radiation inversions form when the ground loses heat via longwave radiation under clear skies at night. This mechanism is fundamentally different from fog-type inversions: the strongest radiation inversions occur on dry, clear nights with light wind, producing large T−Td spreads that cause the dewpoint convergence signal ($S_1$) to return near-zero values. Before this signal was added, the model had a systematic blind spot for radiation inversions — a verified high-intensity calibration event on 18 March 2026 recorded strong odor inside a residence at 6 AM under clear skies with a 9°F dewpoint spread, while the model predicted “Low” (score 9). This event served as a ground-truth benchmark to refine the signal weights ($S_8$) for radiation-driven stability.

The signal employs a three-tier approach with increasing physical fidelity depending on available data. All tiers are active only during the nighttime radiation cooling window (7 PM–8 AM) and weighted by time-of-night (1.0 midnight–6 AM, 0.8 at sunrise transition, 0.7 in evening).

Tier 1: METAR Cloud Ceiling Height + Brutsaert Radiative Cooling (observation periods)

When METAR cloud layer data is available, the system extracts the cloud ceiling height — defined as the lowest broken (BKN, 5–7 oktas) or overcast (OVC, 8 oktas) layer, per aviation convention. This physical measurement distinguishes between low clouds that suppress radiation inversions (warm, strong longwave re-emitters) and high/mid-level clouds that do not (cold, negligible longwave contribution). The previous text-based approach classified METAR “Cloudy” as a single category, missing the critical distinction between low overcast at 1,000 ft (which kills radiation cooling) and altocumulus at 15,000 ft (which barely affects it):

Ceiling HeightRadiation FactorPhysical Basis
< 1,000 ft0.05Stratus/fog deck — warm, nearly blackbody re-emitter
1,000–2,000 ft0.15Low overcast — strong downwelling suppression
2,000–3,000 ft0.25Low-mid cloud deck
3,000–5,000 ft0.45Mid-level — partial suppression
5,000–8,000 ft0.65Mid-level — moderate cooling still occurs
8,000–12,000 ft0.80High-mid (altocumulus) — minor LW impact
> 12,000 ft0.90Cirrus/high cloud — negligible LW contribution

This ceiling factor is combined with the Brutsaert radiative cooling estimate (below) at a 60/40 weighting:

$$S_8 = \left(0.60 \times F_{\text{ceil}} + 0.40 \times R_{\text{signal}}\right) \times w_{\text{night}}$$

Brutsaert (1975) Radiative Cooling Subroutine

All tiers of $S_8$ use the Brutsaert (1975) clear-sky emissivity parameterization to estimate the net longwave radiative cooling rate at the surface. The Brutsaert equation was selected over the older Brunt (1932) empirical formula based on global model intercomparison studies (Morales-Salinas et al. 2023) that found it to be the best-performing emissivity model across all climate zones, with the lowest RMSE and highest coefficient of determination.

The clear-sky atmospheric emissivity is:

$$\varepsilon_{\text{clear}} = 1.24 \left(\frac{e_a}{T_K}\right)^{1/7}$$

where $e_a$ is the vapor pressure at the dewpoint (hPa), computed via August-Roche-Magnus:

$$e_a = 6.112 \exp\!\left(\frac{17.67\, T_{d,C}}{T_{d,C} + 243.5}\right)$$

and $T_K$ is the air temperature in Kelvin.

Two corrections are applied:

1. Nighttime emissivity correction (Li, Jiang & Coimbra 2017): The nocturnal atmosphere is a slightly better longwave emitter than the daytime atmosphere due to the nocturnal temperature profile and moisture redistribution. Li et al. found that nighttime clear-sky emissivity exceeds daytime values by approximately 0.035:

$$\varepsilon_{\text{night}} = \varepsilon_{\text{clear}} + 0.035$$

2. Cloud modification factor (Crawford & Duchon 1999): When METAR cloud fraction ($CF$, 0–1) is available, clouds increase effective emissivity by re-radiating longwave back to the surface. The METAR cloud amounts are converted to fractional coverage (FEW = 0.19, SCT = 0.38, BKN = 0.75, OVC = 1.0) and applied as:

$$\varepsilon_{\text{all}} = \varepsilon_{\text{night}} \times \left(1 + 0.22 \times CF^{2.75}\right)$$

The net radiative cooling rate is then:

$$R_{\text{net}} = \sigma\, T_K^4\, (1 - \varepsilon_{\text{eff}})$$

where $\sigma = 5.67 \times 10^{-8}$ W/m²/K&sup4; is the Stefan-Boltzmann constant. This is normalized to a 0–1 signal with a linear mapping: $R_{\text{net}} < 20$ W/m² → 0 (negligible cooling, e.g., overcast), $R_{\text{net}} > 70$ W/m² → 1.0 (strong cooling, e.g., clear dry night).

The key insight is that the Brutsaert equation reverses the sign of the humidity effect relative to the dewpoint convergence signal ($S_1$): dry air produces lower atmospheric emissivity → higher net cooling → stronger radiation inversion signal. This is physically correct — dry, clear nights are the strongest radiation inversion producers — and eliminates the false-negative blind spot where $S_1$ returned near-zero for exactly the conditions that produce the most dangerous inversions.

Tier 2: Forecast Text + Brutsaert Cooling (forecast-only periods)

When METAR ceiling data is unavailable (future forecast hours), the NWS shortForecast text is parsed into a cloud factor and combined with the Brutsaert cooling signal at 55/45 weighting:

Sky Condition TextCloud Factor
Clear / Fair1.0
Partly Cloudy0.45
Mostly Cloudy / Overcast0.1
Fog / Mist0.3

Tier 3: Forecast Text Only (minimal data)

When neither ceiling height nor temperature/dewpoint data is available, the text-based cloud factor alone is used, scaled by the time-of-night weight.

Signal 9: Inversion Persistence (Weight: 0.20 forecast / 0.15 observation)

Inversions do not break instantaneously when surface conditions change. A radiation inversion that formed under clear skies overnight persists until solar heating generates enough convective mixing to erode the stable layer — typically 1–3 hours after sunrise depending on cloud cover and wind speed. Without this signal, each forecast hour is scored independently: if a METAR observation changes from “Clear” to “Cloudy” at 7 AM (e.g., mid-level altocumulus moving in at sunrise), the sky condition signal drops from 1.0 to 0.1 even though the trapped emissions from overnight are still present at the surface.

The persistence signal examines previous periods (up to 6 hours back) for inversion-forming conditions — clear skies and calm wind during overnight hours (7 PM–7 AM). For each qualifying prior hour, the formation strength is taken as $\min(S_8, S_3)$ (requiring both clear sky and calm wind), then decayed exponentially based on elapsed time:

$$S_9 = \max_{h \in \text{prev}} \left[\min(S_8^{(h)}, S_3^{(h)}) \cdot e^{-0.35\,\Delta t_h}\right]$$

where $\Delta t_h$ is the hours elapsed since period $h$. The decay constant of 0.35 produces a half-life of approximately 2 hours, meaning that a strong inversion signal (1.0) from 4 AM still contributes ~0.25 at 8 AM even if current conditions no longer favor inversion formation. The signal is active only during the morning transition window (4 AM–noon); before 4 AM the current sky condition signal handles detection directly, and after noon solar mixing has had sufficient time to break any residual inversion.

4.4.2 Combined Inversion Score

The system uses dual-mode weighting across nine signals depending on whether real-time observation data is available for the current period. When observation data is present, weights shift to favor measured signals over inferred ones; when unavailable (future forecast hours), the forecast-only weights apply. In both modes, the sky condition ($S_8$) and persistence ($S_9$) signals carry significant weight because they capture radiation inversions that the dewpoint signal misses:

Forecast-only mode (future hours):

$$S_{\text{inv}} = 0.10\,S_1 + 0.10\,S_2 + 0.20\,S_3 + 0.05\,S_4 + 0.15\,S_5 + 0.20\,S_8 + 0.20\,S_9$$

Observation-blended mode (current hour):

$$S_{\text{inv}} = 0.10\,S_1 + 0.10\,S_2 + 0.15\,S_3 + 0.05\,S_4 + 0.10\,S_5 + 0.15\,S_6 + 0.10\,S_7 + 0.10\,S_8 + 0.15\,S_9$$

The rationale for the current weight distribution: the dewpoint convergence signal ($S_1$) was reduced from its original 0.25 to 0.10 because it produces false negatives for radiation inversions (clear, dry nights with large T−Td spreads). The sky condition ($S_8$) and persistence ($S_9$) signals together receive 0.40 in forecast mode, reflecting that radiation inversions are the dominant nocturnal inversion type in the region. In observation mode, visibility ($S_6$) directly measures atmospheric trapping, and $S_8$/$S_9$ are reduced slightly since the measured signals provide more reliable evidence. Both observation-only signals ($S_6$, $S_7$) return 0 when observation data is unavailable, ensuring zero regression in forecast-only mode.

Subject to the wind veto: if wind speed ≥ 8 mph, $S_{\text{inv}} = \min(S_{\text{inv}},\; 0.2)$.

4.4.3 Inversion Amplification Factor

The inversion score is converted to a residual amplification factor that captures concentration enhancement effects not already accounted for by the BVF-dependent dispersion in the transport factor (Section 4.3.3). These residual effects include:

$$I = \begin{cases} 1.0 & \text{if } S_{\text{inv}} < 0.3 \\ 1.0 + (S_{\text{inv}} - 0.3) \times \frac{0.8}{0.7} & \text{if } S_{\text{inv}} \geq 0.3 \end{cases}$$

At maximum inversion confidence ($S_{\text{inv}} = 1.0$), the residual amplification factor reaches approximately 1.8×. This was increased from a previous value of 1.35× following the 18 March 2026 calibration event, which indicated that radiation inversions — where emissions accumulate under a clear-sky stable layer for 6–10 hours overnight — produce substantially higher ground-level concentrations than the earlier factor captured. The combined effect of both the inversion factor and the stability-dependent BVF dispersion parameterization in the transport factor at $S_{\text{inv}} = 1.0$ produces a net concentration increase of approximately $1.8 \times 1.8 \approx 3.2\times$ relative to neutral conditions at moderate distance, which better matches the magnitude of observed odor intensification during strong nocturnal inversions.

4.5 Diurnal Factor (D)

The diurnal factor accounts for the daily cycle of atmospheric mixing driven by solar heating. Traditional dispersion models apply strong daytime suppression (mixing height increases → dilution) because they focus on elevated point-source emissions that are lofted and dispersed by convective mixing. However, the open-air area sources at these facilities — treated wood storage yards spread across acres of ground-level surface — are the primary daytime odor driver. These area sources have zero effective stack height: emissions originate at the surface where people breathe, and even vigorous convective mixing does not eliminate ground-level impact the way it disperses emissions from elevated stacks. Simultaneously, solar heating maximizes area source off-gassing precisely when convective mixing is strongest.

The diurnal factor uses solar-relative timing rather than fixed clock hours. Phase boundaries are defined relative to the astronomically computed sunrise ($t_r$), sunset ($t_s$), and solar noon ($t_n$) for the actual forecast date, so summer evenings (sun still strong at 7 PM CDT) correctly maintain higher scores, and winter mornings (dark until nearly 7 AM) do not apply solar-driven factors prematurely:

PhaseTime WindowFactorPhysical Basis
Sunrise$t_r - 1$ to $t_r + 1$1.0Boundary layer collapse + area surfaces beginning to warm
Late morning$t_r + 1$ to $t_n - 1$0.5Convective mixing ramps up; solar heating of surfaces increasing
Midday$t_n - 1$ to $t_n + 1$0.55Peak solar heating of area sources partially offsets peak mixing
Afternoon$t_n + 1$ to $t_s - 1$0.65Peak surface temps; convective mixing beginning to decline
Evening$t_s - 1$ to $t_s + 3$0.8Surfaces still warm; vertical mixing drops rapidly post-sunset
Night$t_s + 3$ to $t_r - 1$0.6Surfaces cool; inversions may trap residual and kiln emissions

The early morning maximum reflects the well-documented phenomenon of boundary layer collapse at sunrise, combined with area sources that are beginning to absorb solar radiation. The midday factor (0.55) is notably higher than the 0.3 used in traditional point-source models, reflecting the physical reality that ground-level area sources maintain significant impact even during peak convective mixing. The afternoon factor (0.65) captures the critical period when surface temperatures peak while atmospheric mixing begins to subside.

Using solar-relative phases ensures the model adapts correctly across the year. For example, at Alexandria’s latitude the “Evening” phase (factor 0.8) extends to approximately 11:13 PM on the summer solstice (sunset ~8:13 PM + 3 hours) but only to ~8:07 PM on the winter solstice (sunset ~5:07 PM + 3 hours). This ~4-hour seasonal swing in daylight duration (10h in winter vs 14h in summer) produces meaningful differences in when thermal retention and source weighting transitions occur, and the model captures this without manual seasonal adjustments.

4.6 Humidity Factor (H)

Humidity affects odor perception and transport through multiple mechanisms, including moisture absorption of water-soluble compounds and enhanced olfactory sensitivity in humid air:

Relative HumidityFactorMechanism
≥ 97%1.30Fog/heavy mist — odor compounds adsorb to droplets
93–97%1.15Near-fog; enhanced dissolution of soluble compounds
80–93%1.00Baseline (typical Louisiana humidity)
50–80%0.90Below-average humidity reduces transport
< 50%0.75Unusually dry; rapid evaporative dispersion

The baseline of 1.0 at 80–93% humidity is calibrated to Louisiana’s climate, where high humidity is the norm rather than an exceptional condition. An additional 10% multiplier is applied to Pineville’s kiln emissions at humidity ≥ 93%, reflecting the high water solubility of acetaldehyde, formaldehyde, and methanol.

4.7 Terrain Factor (G)

The terrain factor models the influence of topography on pollutant transport. Unlike flat-terrain dispersion models that treat the landscape as a uniform surface, this system explicitly accounts for the geological features of the Alexandria–Pineville area that profoundly affect where emissions accumulate and how they move.

4.7.1 Geological Context

The study area is defined by the Red River valley, a north-south alluvial floodplain carved through central Louisiana. The valley floor sits at 70–78 ft elevation and runs roughly north-south, creating a natural channel for airflow. On the east side of the river, the Pineville bluff rises sharply to 100–105 ft — a 30-foot escarpment that has direct consequences for pollutant behavior:

4.7.2 Elevation Model

The system uses inverse-distance-weighted (IDW) interpolation from 50+ elevation control points queried from the USGS National Map Elevation Point Query Service (1/3 arc-second DEM) to estimate terrain height at any latitude-longitude coordinate:

$$z(\mathbf{x}) = \frac{\sum_{i=1}^{25} w_i\, z_i}{\sum_{i=1}^{25} w_i}, \quad w_i = \frac{1}{d(\mathbf{x}, \mathbf{x}_i)^2}$$

where $z_i$ is the known elevation at control point $\mathbf{x}_i$ and $d(\mathbf{x}, \mathbf{x}_i)$ is the Haversine distance. The control points are concentrated along the key terrain features:

FeatureElevationControl PointsSignificance
Valley floor (Alexandria)72–78 ft10Natural pooling basin for dense emissions
Pineville bluff95–105 ft7Sharp escarpment, drainage source
Western rise80–82 ft4Gradual slope, partial containment
East of Pineville75–95 ft4Transition zone back to lower terrain

The IDW approach with squared distance weighting produces a smooth elevation surface that captures the dominant features — particularly the sharp Pineville bluff — while remaining computationally inexpensive (a single pass through 25 points per receptor location).

4.7.3 Drainage Flow Modeling

Under calm or inversion conditions (wind ≤ 5 mph, inversion score > 0.4), the vertical atmospheric mixing that normally disperses pollutants upward is suppressed. In this regime, gravity becomes the dominant transport mechanism for heavy compounds. Dense vapors behave like cold air in a katabatic flow — sliding downhill along the terrain surface and pooling in topographic depressions.

For receptor points below the plant elevation ($\Delta z < -10$ ft):

$$G_{\text{drain}} = 1.0 + \min\!\left(0.4,\; \frac{|\Delta z|}{75}\right)$$

This produces a maximum 40% concentration enhancement for points deep in the valley relative to the source. For the Pineville facility specifically, drainage from the bluff receives a 1.5× enhancement when the elevation difference exceeds 20 ft, reflecting two compounding physical mechanisms:

  1. Gravitational drainage: The 30-ft bluff creates a steep gradient that accelerates dense vapor flow downhill. Unlike gentle slopes where drainage is sluggish, the Pineville escarpment produces a concentrated drainage stream that carries heavy compounds (diesel, cresol) efficiently into the valley.
  2. Valley pooling: The valley floor acts as a topographic bowl. Once dense emissions reach the bottom, they have nowhere to go — the western rise (80–82 ft) and the river channel create partial containment. During inversions, this pooled mass cannot disperse vertically, leading to sustained high concentrations at breathing height.

For receptor points above the plant elevation ($\Delta z > 15$ ft), the terrain factor is reduced to 0.5, reflecting the shielding effect of higher ground. Points slightly above ($\Delta z$ 5–15 ft) receive a modest reduction to 0.8.

4.7.4 Valley Channeling

Under light-to-moderate wind conditions (≤ 10 mph), winds aligned with the north-south Red River valley axis (0° or 180° ± 30°) receive a channeling boost of 1.15× for valley-floor receptor points. This models the well-documented tendency of valley topography to constrain and accelerate airflow along the valley axis, concentrating pollutants into a narrower corridor rather than allowing lateral dispersion across the full wind rose.

4.7.5 Fenceline Exclusion

The model excludes a 150-meter radius around each plant center from community odor scoring. This fenceline exclusion serves two purposes: (1) the model is designed to predict community exposure, not occupational exposure within the facility boundary, and (2) the near-source concentration field is dominated by turbulent wake effects from buildings and process equipment that the simplified dispersion model cannot capture. Scores within the exclusion zone return zero, and the user interface displays “Facility Zone — Not Tracked” when a user places a pin inside the fenceline.

4.7.6 Terrain Data Source and Resolution

Elevation values are queried from the USGS National Map Elevation Point Query Service (EPQS), which provides access to the 1/3 arc-second (~10 m) Digital Elevation Model. The current grid uses 50+ control points with increased density in critical areas:

4.7.7 Current Limitations of the Terrain Model

The IDW elevation model captures the dominant geological features but does not account for:


5. Chemical-Specific Plume Geometry

5.1 Compound-Specific Transport Characteristics

Beyond the scalar intensity calculation, the system generates chemical-specific plume geometries for map visualization. Each compound class has distinct transport characteristics based on molecular weight, volatility, and release height:

Table 6: Chemical Plume Parameters

Compound Classσ MultiplierDistance MultiplierPhysical Basis
Naphthalene1.30.7MW 128; pools in calm air, wide/short plume
Cresol1.20.65MW 108; dense, strong odor at low concentration
Diesel1.40.5MW ~170; heaviest, widest spread, shortest range
Kiln gases0.71.4MW 30–44; buoyant, narrow/far, elevated release

These multipliers are applied to the base plume geometry parameters, producing visually distinct plume shapes on the map that communicate the different transport behaviors of each compound class.

5.2 Chemical-Specific Heatmap Rendering

Each chemical overlay is rendered as a canvas heatmap using the compound’s designated color. The per-pixel score is computed using the standard prediction model, with the compound’s transport parameters (Table 6) limiting the effective range. The resulting score is then modulated by the ground-level fraction $F_g$ (Section 5.4), so that elevated stack releases appear faint when emissions stay aloft and progressively opaque during inversions and fog.

This produces visualizations that directly communicate compound-specific transport behavior: diesel appears as a dense, concentrated cloud near the source, while kiln gases appear as a diffuse, far-reaching haze.

5.3 Plume Polygon Generation

For wind-driven conditions, plume polygons are generated as fan-shaped sectors:

  1. The downwind direction is computed as $(D_{\text{wind}} + 180°) \bmod 360°$
  2. Angular rays are cast from $-3\sigma$ to $+3\sigma$ in 5° increments
  3. Each ray extends to a distance weighted by the Gaussian angular factor
  4. Destination points are calculated using the Haversine destination formula
  5. The polygon is closed through the plant location

Under calm conditions, a circular polygon is generated with a radius that extends under inversion conditions:

$$r = r_{\text{base}} \times (1 + S_{\text{inv}} \times k)$$

where $k = 0.5$ for ground-level compounds and $k = 0.15$ for stack-released kiln gases (which are less affected by surface-level inversions).

5.4 Ground-Level Fraction

For elevated stack releases (e.g., kiln gases released at 30–60 m height), only a fraction of the plume reaches breathing height. The ground-level fraction $F_g$ estimates the proportion of emissions that affect ground-level receptors, modulating both the chemical overlay opacity and the effective odor contribution. Ground-level releases always have $F_g = 1.0$.

For elevated releases, the base fraction starts at $F_g = 0.35$ (clear-sky, moderate-wind conditions) and is modified by several atmospheric factors:

Table 8: Ground-Level Fraction Modifiers for Elevated Releases

FactorConditionAdjustmentPhysical Basis
Inversion trapping$S_{\text{inv}} > 0.3$$+(S_{\text{inv}} - 0.3) \times 0.75$ Inversion lid traps elevated emissions, forcing plume downward; strong inversions can add up to +0.525
Wind-induced mixing$v_w \ge 3$ mph$+0.05$ to $+0.15$ Mechanical turbulence from moderate winds enhances vertical mixing toward ground
Molecular weightMW > 44$+$up to $0.15$ Heavier-than-air compounds settle toward ground due to density-driven sinking velocity
Humidity/fogHumidity > 70%$+0.05$ to $+0.20$ Fog droplets absorb water-soluble gases; wet deposition pulls plume to surface
Temperature stabilityTemp < 50°F / > 85°F$\pm 0.08$ Cool air suppresses buoyancy (plume sinks); hot air enhances buoyant rise (plume stays aloft)

The final fraction is clamped to $[0.1, 1.0]$. Under worst-case conditions (strong inversion, high humidity, heavy compound, moderate wind), elevated releases can approach $F_g \approx 1.0$, meaning nearly all emissions reach breathing height—reflecting the “fumigation” events that residents experience during temperature inversions.


6. Visualization and User Interface

6.1 Interactive Map

The system renders predictions on an interactive Leaflet.js map with OpenStreetMap base tiles. Key features include:

6.2 Heatmap Plume Visualization

The system renders all plume overlays — both the combined odor plume and individual chemical overlays — using a canvas-based heatmap approach. This provides spatially accurate, per-pixel intensity rendering across the entire map.

6.2.1 Combined Odor Heatmap

The combined heatmap renders a per-pixel combined score from both facilities at every point on the map. A canvas tile overlay samples a grid of points (8-pixel resolution per cell, with GPU-accelerated CSS blur for smoothing), computes computeCombinedScore at each point — which takes the maximum contribution from either plant accounting for emission, transport, inversion, time-of-day, humidity, and terrain factors — and colors each cell according to the score.

Why the heatmap matters: The heatmap’s key advantage is that it shows interaction zones where emissions from both plants overlap. Because both facilities are within 4 km of each other and often share the same downwind corridor, their plumes frequently converge over residential areas. The heatmap computes the actual combined impact at each point, identifying zones of modeled concentration overlap that per-plant visualizations cannot show.

Color is determined by smooth interpolation between risk level colors using a square-root curve that accelerates transitions near level boundaries:

$$\text{color}(s) = \text{lerp}\!\left(\text{color}_{\text{prev}},\; \text{color}_{\text{next}},\; \sqrt{\frac{s - s_{\text{prev}}}{s_{\text{next}} - s_{\text{prev}}}}\right)$$

This eliminates the hard color boundaries that would otherwise appear at score thresholds (e.g., the jump from green to orange at score 15), producing a continuous color field that better reflects the gradual nature of real-world odor intensity transitions.

Opacity uses a two-phase function to keep the Low-intensity zone (green) visually subtle while ensuring Moderate-and-above zones (orange, red) are immediately prominent:

$$\alpha(s) = \begin{cases} 0.22 + \frac{s}{15} \times 0.12 & \text{if } s \leq 15 \\ 0.40 + \frac{s - 15}{85} \times 0.50 & \text{if } s > 15 \end{cases}$$

This design choice reflects community feedback: green zones should be visible but not alarming (scores 0–15 represent conditions where odor is unlikely to be noticeable), while orange and red zones should demand attention because they represent conditions where residents may be affected.

6.2.2 Chemical-Specific Heatmaps

Individual chemical overlays use the same canvas heatmap rendering approach, but colored with each compound’s distinct color and scored using compound-specific transport parameters (sigma multiplier, distance multiplier, and release height). Each chemical heatmap computes per-pixel intensity modulated by the ground-level fraction $F_g$ (Section 5.4), so elevated stack releases appear faint when emissions stay aloft and progressively opaque during inversions and fog.

Enabling any chemical overlay automatically disables the combined plume to avoid visual clutter and allow clear inspection of individual compound dispersion patterns.

6.3 Ground-Level Fraction in Visualization

Chemical heatmap intensity is modulated by the ground-level fraction $F_g$ (Section 5.4). This means that elevated stack releases appear faint on the map during clear, windy conditions (when most emissions stay aloft) and progressively opaque during inversions and fog (when emissions are trapped at breathing height). The visualization thus reflects not just where chemicals disperse, but what people actually smell at ground level.

6.4 Scoring Visualization

A real-time SVG arc gauge displays the current odor score (0–100) with animated transitions, accompanied by current conditions (temperature, humidity, wind, inversion status) and a text description of the forecast.

6.5 Odor Level Classification

Score RangeLevelColorDescription
0–15LowGreen (#4CAF50)Minimal odor expected
16–40ModerateOrange (#FF9800)Noticeable odor possible
41–65ElevatedRed (#F44336)Significant odor likely
66–100HighDark Red (#B71C1C)Strong odor expected

7. Field Calibration and Validation

7.1 Calibration Methodology

The system was calibrated through two complementary approaches: (1) iterative comparison of model predictions against real-time odor observations from community members in the Alexandria–Pineville area, and (2) remote sensing validation using publicly available aerial imagery (Google Maps/Earth) and Street View photography to identify and characterize emission sources not fully captured by permit data alone.

7.2 Key Calibration Adjustments

Naphthalene floor correction (March 2026): Initial model versions used a vapor pressure floor of 0.2 at low temperatures, which underestimated odor at 40–45°F. Field observations during a 42°F morning with confirmed odor detection led to raising the floor to 0.35, reflecting the reality that creosote-soaked surfaces (railroad ties, treatment equipment) maintain significant naphthalene evaporation even at low temperatures due to the large surface area of saturated material.

Distance decay exponent: The decay exponent was adjusted from an initial value of 2.0 (inverse-square law) to 1.8 based on observations of detectable odor at distances greater than predicted by the theoretical model. This adjustment accounts for the fact that complex terrain and building-wake effects can maintain higher concentrations at distance than idealized open-terrain models predict.

Open-air storage area source identification (March 2026): Remote sensing validation of publicly available aerial imagery indicated that both facilities maintain extensive open-air storage of treated wood products that constitutes a dominant emission source not adequately captured by the original fugitive emissions model. At Alexandria, acres of dark creosote-treated railroad ties are visible in satellite imagery, along with untreated (lighter-colored) tie stacks. At Pineville, dozens of treated telephone pole bundles are visible as dark stacks across several acres, alongside the kiln structure with rooftop exhaust vents. This discovery led to the multi-source emission model with solar surface temperature calculations and diurnally shifting source weights, substantially improving daytime prediction accuracy.

Diurnal factor recalibration: The original model applied a 0.3 suppression factor during midday hours, consistent with traditional point-source dispersion modeling where convective mixing disperses emissions vertically. After incorporating the area source model, the midday factor was raised to 0.55 to reflect that ground-level area sources (zero stack height) maintain significant impact even during peak convective mixing, and that solar heating maximizes area source emissions precisely when atmospheric mixing is strongest.

Radiation inversion calibration (March 2026): A verified calibration event on 18 March 2026 — odor detectable inside a residence at 6 AM while the model predicted “Low” (score 9) — exposed a systematic blind spot: the inversion detector relied on dewpoint convergence (T−Td < 2°F) as its primary signal, but the event occurred under clear skies with a 9°F dewpoint spread, conditions that produce strong radiation inversions through longwave surface cooling rather than moisture-driven stability. Two new signals were added: (1) a sky condition signal ($S_8$) that detects clear-sky radiative cooling at night using NWS shortForecast text, and (2) an inversion persistence signal ($S_9$) that carries overnight inversion-forming conditions forward through the morning transition with exponential decay (half-life ~2 hours), preventing the score from collapsing when a METAR observation changes from “Clear” to “Cloudy” while the trapped layer is still present. The inversion amplification factor was also increased from 1.35× to 1.8× to reflect the prolonged emission accumulation under nocturnal radiation inversions. Signal weights were redistributed to reduce the dewpoint signal from 0.25 to 0.10 and allocate 0.40 total weight to the new sky condition and persistence signals in forecast mode. Subsequently, $S_8$ was upgraded from text-only classification to a three-tier physics-based approach: (1) METAR cloud ceiling height extraction distinguishes low overcast (which suppresses radiation cooling) from high/mid-level cloud (which does not), (2) the Brutsaert (1975) atmospheric emissivity equation estimates net longwave cooling rate from temperature and dewpoint, with a nighttime correction from Li, Jiang & Coimbra (2017) and cloud modification from Crawford & Duchon (1999), and (3) the ceiling measurement and Brutsaert cooling signal are combined at 60/40 weighting to produce a physically grounded radiation inversion probability rather than relying on text parsing of ambiguous METAR descriptions.

Solar position model (March 2026): All daytime/nighttime determinations were replaced with astronomically computed sunrise and sunset times using the Spencer (1971) solar declination formula. The solar heating boost, source weighting shifts, and diurnal factor phases are now defined relative to actual solar position rather than fixed clock hours. Additionally, the solar heating boost is scaled by the normalized solar elevation angle (0 at horizon, 1 at solar noon), so surface temperature estimates ramp gradually after sunrise rather than switching on abruptly. This corrects a ~4-hour seasonal swing in daylight (from ~10h at winter solstice to ~14h at summer solstice) that the fixed-hour model could not capture.

7.3 Validation Results

The system’s predictive output has shown high correlation with community-reported odor events, confirmed through real-time observations. The multiplicative factor architecture has proven particularly effective at capturing the interaction effects between meteorological variables — for example, correctly predicting elevated odor during cool, calm, humid mornings with temperature inversions, even when individual factors (e.g., low temperature reducing emissions) might suggest otherwise. The area source model improved prediction accuracy during hot afternoon periods, which the original model had systematically underestimated due to excessive daytime suppression that did not account for solar-driven ground-level off-gassing.


8. Technical Architecture

8.1 Application Stack

The application is implemented as a static web application with no server-side processing requirements:

LayerTechnologyPurpose
UI FrameworkPreact + htmLightweight reactive components
MappingLeaflet.jsInteractive map rendering
Data SourcesNWS API + Open-Meteo APINWS: KAEX observations, hourly forecasts, AFD; Open-Meteo: ECMWF/GFS model forecasts for multi-model blending
Data IngestionCloudflare Cron Trigger (every 30 min)Sole external API consumer — fetches 4 sources in parallel, archives all data to D1
Data StorageCloudflare D1 (SQLite)NWS forecasts, Open-Meteo forecasts, observations, AFD text, odor reports; server-side blending at serve time
API LayerCloudflare WorkerServes /api/* endpoints from D1 — no NWS proxy
DeploymentStatic hostingZero-infrastructure operation

8.2 Performance Considerations

The system achieves responsive performance through:


9. Discussion

9.1 Significance of AI-Assisted Development

This project demonstrates a paradigm shift in environmental monitoring tool development. Traditional atmospheric dispersion modeling requires:

By contrast, this system was developed through iterative dialogue with an AI assistant, translating domain knowledge from environmental permits, atmospheric science textbooks, and field observations into working code. The AI synthesized principles from:

This cross-disciplinary synthesis, performed through natural language interaction, enabled a non-specialist to develop a tool with observable correlation to field conditions.

9.2 Model Classification and Scope

To preempt misinterpretation of this system’s claims, it is important to situate it within the established hierarchy of atmospheric dispersion models:

  1. Teaching tools (e.g., NOAA READY Gaussian Plume): Single-equation steady-state Gaussian models with no plume rise, no meander, no terrain effects. NOAA’s own documentation describes READY as “a teaching tool using a very simple model” and recommends HYSPLIT for all serious dispersion work.
  2. Screening and community tools (this system): Incorporates Briggs plume rise, AERMOD meander transition, BVF stability-dependent dispersion, log-law wind profile, terrain-aware transport, and chemical-specific volatilization physics. Produces relative odor likelihood scores, not regulatory-grade concentrations. Uses real-time NWS observations and forecasts. Validated against field odor reports.
  3. Regulatory models (AERMOD, CALPUFF, HYSPLIT): Full implementation of boundary-layer turbulence parameterization, building downwash, hour-by-hour meteorological processing, source-specific emission inventories, and concentration output in µg/m³. Require professional meteorological input (AERMET-processed surface and upper-air data), professional consulting for configuration and interpretation, and expert operation. Used for permit applications, NAAQS compliance demonstrations, and enforcement actions.

This system does not claim to replace regulatory modeling. It occupies a niche that regulatory models do not serve: real-time, publicly accessible, condition-specific odor forecasting for community use. No regulatory model provides this capability because they are designed for retrospective compliance analysis, not prospective community alerts. The relevant question is not “is this as good as AERMOD?” but “does this provide useful information that communities would not otherwise have?”

9.3 Limitations

9.3.1 Fundamental Limitations of the Gaussian Approach

  1. Odor intermittency and the peak-to-mean problem: This is the most significant theoretical limitation of any Gaussian plume application to odor prediction. Gaussian models compute time-averaged concentration fields—smooth, bell-shaped distributions that represent the statistical mean over the averaging period (typically one hour). Human odor perception, however, responds to instantaneous concentration fluctuations: brief, intense “whiffs” produced by turbulent eddies that transport high-concentration parcels intermittently (Barclay & Scire, 2021). The ratio of instantaneous peak concentration to the time-averaged mean (the peak-to-mean ratio) can range from 5 to 50 or more near the source and at plume edges, depending on turbulence intensity, averaging time, and distance. This means a Gaussian model predicting a “low” average concentration may coexist with intermittent peaks well above the odor detection threshold.

    Implication for this system: The model will tend to underpredict odor detection at plume edges (where intermittent whiffs are detectable but the time-averaged concentration is low) and overpredict the spatial uniformity of odor within the plume (the real experience is patchy, not smooth). This is partially mitigated by two factors: (a) the system’s odor scores represent likelihood, not concentration, so a moderate score in a region with high intermittency still communicates meaningful odor risk; and (b) the hourly resolution naturally captures the conditions (wind speed, stability) that correlate with high peak-to-mean ratios (stable, low-wind conditions produce both higher average concentrations and higher intermittency, so the score rises when whiffs are most likely). A formal peak-to-mean correction factor, as implemented in fluctuating plume models, remains a candidate for future refinement.
  2. Steady-state assumption: Each hourly prediction assumes meteorological conditions are constant throughout the hour. Rapid wind shifts, frontal passages, or convective onset can produce transient odor events not captured by the hourly resolution. The system cannot model the residual plume from the previous hour or the transition between two different wind regimes within a single hour.
  3. No building downwash or obstacle effects: The model does not account for aerodynamic wakes behind buildings, tree lines, or other obstacles. In the immediate vicinity of the facilities, building-induced turbulence can trap emissions in recirculation zones or loft them above the modeled plume trajectory. AERMOD’s BPIP preprocessor handles this; this system does not.
  4. No chemical transformation during transport: The model treats all emitted compounds as chemically inert during transport. In reality, some creosote VOCs undergo photolysis and hydroxyl radical reaction, particularly during daytime. This means daytime predictions may slightly overestimate concentrations at long range—a conservative bias.

9.3.2 Data and Calibration Limitations

  1. Unknown emission rates: The system does not know actual emission rates ($Q$) from either facility. LDEQ permits specify allowable limits, not real-time emissions, and actual emissions vary with production rate, feedstock, equipment condition, and operator behavior. Without $Q$, the model cannot compute absolute concentrations and instead produces relative scores. This means the system can identify when and where odor is most likely but cannot quantify how much pollutant is present in regulatory units (µg/m³). Any claim of “prediction” should be understood as “relative likelihood given atmospheric conditions,” not “measured pollutant concentration.”
  2. Single-point meteorological data: The system uses a single NWS grid point and KAEX ASOS station, which may not capture microscale variations in wind speed and direction across the study area. Wind channeling through the Red River valley, local thermal circulations, and terrain-induced eddies can produce conditions that differ from the single measurement point.
  3. Forecast-dependent accuracy (future hours): Predictions beyond the current hour are only as accurate as the underlying NWS forecast. Past and current hours now use KAEX observations, and the D1 archive enables systematic measurement of forecast error to inform future bias corrections.
  4. Qualitative validation: Field validation relies on subjective odor perception rather than instrumental VOC measurements. Odor detection thresholds vary between individuals by up to an order of magnitude. The system has not been validated against collocated chemical sensors. Until community odor reports are systematically collected and correlated with model predictions (a planned capability), the validation remains anecdotal rather than statistical.

9.3.3 What This System Is and Is Not

To be explicit about scope: this system is a weather-driven odor likelihood heuristic, not an atmospheric concentration model. It answers the question “given current and forecasted weather conditions, how likely is odor to be detectable in my area?” It does not answer “what is the ambient concentration of naphthalene at this receptor?” These are fundamentally different questions, and conflating them would be misleading.

The system is:

The system is not:

Treating the output as more than it claims to be would be a misuse. Treating it as less than it provides—dismissing weather-driven odor likelihood as mere “vibes”—ignores that the atmospheric conditions governing odor transport (wind speed, direction, stability, humidity, inversion state) are measurable physical quantities with well-understood effects on pollutant dispersion, and that this system integrates those quantities through equations with documented provenance (Appendix C).

9.4 Community Empowerment Through Condition-Specific Visualization

A central insight of this project is that visualized, condition-specific data transforms community advocacy. The statement “it stinks” is a subjective complaint that is easy to dismiss. The statement “a temperature inversion trapped ground-level creosote emissions under 85% humidity with 3 mph winds from the south, producing elevated odor conditions at our school during morning drop-off” is an evidence-based observation that demands a response.

The system’s heatmap plume visualization makes this transformation tangible. When a resident sees a dense red heatmap engulfing their neighborhood during a morning inversion, they are not merely experiencing a bad smell—they are seeing the atmospheric physics that explain why it smells, how far the odor extends, and which conditions are responsible. This shifts the conversation from subjective perception to objective meteorology.

This data-driven approach empowers communities in several ways:

In this sense, the system serves not as an accusation but as a shared language—one grounded in atmospheric science rather than subjective perception—through which communities and industry can work together to minimize odor impact by accounting for the weather realities that neither party controls.

9.5 Future Directions

Recently implemented:

Planned:


10. Conclusion

This work presents a novel, accessible approach to community-level industrial odor forecasting that combines real-time weather data, chemical-specific volatilization physics, and terrain-aware atmospheric dispersion modeling into a unified, browser-based prediction system. The heatmap plume visualization—with per-pixel scoring, chemical-specific transport parameters, and ground-level fraction modulation—communicates not just where odor may occur, but why and how intensely, grounding subjective experience in atmospheric physics.

The system’s strong correlation with community-reported odor events validates the underlying atmospheric model as a tool for identifying conditions that amplify localized odor perception. This correlation demonstrates the efficacy of using publicly available atmospheric and permit data as a proxy for neighborhood-level impact. By making condition-specific data accessible to non-specialists, the system supports evidence-based community engagement with environmental conditions. When residents can articulate the meteorological conditions behind their experience, and when industry can anticipate those conditions proactively, the path opens from confrontation to collaborative stewardship.

The project demonstrates that large language models can serve as effective collaborators in developing sophisticated environmental monitoring tools, synthesizing domain knowledge across multiple scientific disciplines and translating it into functional software without requiring specialized expertise or expensive infrastructure. This broader accessibility of environmental science tools has significant implications for community environmental monitoring and the goal of informed coexistence between industry and the communities they operate within.


Appendix A: Mathematical Notation Summary

SymbolDefinitionUnits
$E$Emission factordimensionless
$T$Transport factordimensionless
$I$Inversion amplification factordimensionless
$D$Diurnal (time-of-day) factordimensionless
$H$Humidity factordimensionless
$G$Terrain (geographic) factordimensionless
$T_F$Ambient air temperature°F
$T_{\text{eff}}$Effective surface temperature (solar-heated)°F
$\Delta T_{\text{solar}}$Solar heating boost (peak, at solar noon)°F
$S_I$Normalized solar intensity (elevation angle / peak elevation)0–1
$\alpha$Solar elevation angledegrees
$\delta$Solar declinationradians
$\omega_0$Hour angle at sunrise/sunsetdegrees
$t_r, t_s, t_n$Sunrise, sunset, solar noon (local time)hours
$v_w$Wind speedmph
$\sigma$Gaussian plume spread parameterdegrees
$d$Haversine distancekm
$d_{\max}$Maximum transport distancekm
$s$Atmospheric stability parameter0–1
$\beta(s)$Stability-dependent distance decay exponentdimensionless
$M_s$Stability-dependent range multiplierdimensionless
$N$Brunt-Väisälä frequencys$^{-1}$
$\sigma_z$Vertical dispersion parameterm
$\sigma_w$Vertical velocity variancem/s
$\eta(s)$Stability-dependent dilution efficiencydimensionless
$S_6$Visibility signal (observation-only)dimensionless
$S_7$Pressure trend signal (observation-only)dimensionless
$P$Barometric pressuremb (hPa)
$V$Visibilitymiles
$\Delta\theta$Angular offset from downwinddegrees
$S_{\text{inv}}$Inversion confidence score0–1
$\Delta z$Elevation difference (receptor – source)ft
$B_I$Base intensitydimensionless
$F_g$Ground-level fraction (proportion of elevated emissions reaching breathing height)0–1

Appendix B: Elevation Control Points

The terrain model uses 50+ elevation control points queried from the USGS National Map Elevation Point Query Service (1/3 arc-second DEM, ~10 m resolution), spanning the Alexandria–Pineville study area. Elevations range from 65 ft (creek drainages southeast of Pineville) to 113 ft (Pineville bluff ridge). Point density is highest along the bluff edge and in the Pineville HS drainage basin, where terrain-driven pooling significantly affects ground-level concentrations. Full coordinates and elevations are available in the application source code.

Appendix C: Model Validation Against Published Literature

This appendix subjects the system’s atmospheric dispersion mathematics to critical review against established literature, identifying every deliberate departure from textbook Gaussian plume theory and evaluating whether each simplification is defensible. The goal is full transparency: a reader should be able to identify exactly where this model departs from reference formulations, understand why it departs, and judge for themselves whether the trade-offs are acceptable.

C.1 Reference Formulation

The standard Gaussian plume equation for a continuous, elevated point source, as formulated by Pasquill (1961) and codified in Turner (1970; revised 1994), takes the “complete” form:

$$C = \frac{Q}{u} \cdot \frac{f}{\sigma_y \sqrt{2\pi}} \cdot \frac{g_1 + g_2 + g_3}{\sigma_z \sqrt{2\pi}}$$

where the terms are:

At ground level ($z = 0$), the direct and reflection terms collapse: $g_1 = g_2 = \exp[-H^2/(2\sigma_z^2)]$, yielding the simplified ground-level centerline equation:

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

The factor of 2 in the numerator (from $g_1 + g_2 = 2g_1$ at $z = 0$) cancels with the denominator’s $2\pi$, producing $\pi$ in the denominator. This is the “mirror-image doubling” that any ground-level receptor inherits from the reflection boundary condition.

All subsequent comparisons reference this formulation. Section numbers refer to the main body of this document.

C.2 Departure Analysis

C.2.1 Relative Scoring vs. Absolute Concentration

Departure: The standard equation computes concentration $C$ in µg/m³ from a known emission rate $Q$ (g/s). This system does not compute absolute concentration. Instead, it produces a dimensionless odor score (0–10) by combining multiplicative atmospheric factors (Section 4.1). The $Q / (u\,\sigma_y\,\sigma_z)$ prefactor is absent.

Criticism: Without absolute concentration, the model cannot be compared against regulatory ambient air quality standards (e.g., NAAQS), cannot produce isopleth maps in standard units, and cannot be validated against instrumental VOC measurements. The dimensionless score has no physical units and its mapping to human odor perception depends on calibration that may not generalize.

Justification: The emission rate $Q$ for these facilities is unknown. LDEQ permits specify allowable emission limits, not actual emission rates, and actual emissions vary with production rate, feedstock moisture, kiln temperature, and equipment condition. Fabricating a $Q$ value would introduce false precision. The relative scoring approach acknowledges this limitation honestly: the model predicts when and where odor is most likely, not how much pollutant is present. For community odor forecasting—where the question is “will I smell it tonight?” rather than “is the ambient concentration above 35 ppb?”—a well-calibrated relative score is more useful than an absolute concentration derived from an uncertain $Q$. The system’s field validation (Section 7) demonstrates that the relative approach correctly identifies real odor events.

C.2.2 Ground-Reflection Factor ($g_2$ Term)

Departure: The ground-level fraction calculation (Section 5.4) uses $F_g = \exp[-H_e^2 / (2\sigma_z^2)]$, which corresponds to the $g_1$ term alone. The $g_2$ term, which equals $g_1$ at $z = 0$ and doubles the ground-level concentration, is not explicitly included. The model therefore underestimates the ground-contact fraction by a factor of approximately 2 relative to the textbook equation.

Criticism: The factor-of-2 omission is a systematic bias. In a model claiming to be “AERMOD-inspired,” omitting the fundamental ground-reflection boundary condition is a significant departure from established dispersion physics.

Justification: The factor of 2 is a multiplicative constant that applies uniformly across all receptor locations and atmospheric conditions. In a relative scoring system without absolute concentration, any constant multiplicative factor is absorbed into the base intensity calibration ($B_I$ in Section 4.2). Including the factor of 2 would simply halve the calibrated base intensity, producing identical final scores. The ground-level fraction $F_g$ is used as a modulator (0–1) on the chemical overlay, not as an absolute concentration. As long as the ratio of $F_g$ values between different atmospheric conditions is preserved—which it is, since the factor of 2 cancels in any ratio—the relative prediction is unaffected.

C.2.3 Inversion-Lid Reflections ($g_3$ Term)

Departure: The textbook $g_3$ term is an infinite series of Gaussian reflection terms between the ground and an inversion lid at height $L$. This model replaces the series with an empirical inversion boost: when the inversion confidence score $S_{\text{inv}} > 0.3$, the ground-level fraction receives an additive increase of $(S_{\text{inv}} - 0.3) \times 0.4$, up to approximately +0.28 (Section 5.4).

Criticism: The empirical boost has no direct physical relationship to the inversion lid height $L$. It cannot reproduce the distance-dependent behavior of the $g_3$ series, which becomes more important as $\sigma_z$ approaches $L$ (i.e., at greater downwind distances or under stronger instability). The linear scaling with $S_{\text{inv}}$ is a crude approximation of a nonlinear physical process.

Justification: The $g_3$ series requires the inversion lid height $L$ as input. This system has no direct measurement of $L$: the inversion score (Section 4.4) indicates the presence and approximate strength of a temperature inversion, but not its altitude. Without $L$, the $g_3$ series cannot be evaluated. The empirical approach captures the first-order physical effect—inversions trap emissions and increase ground-level concentration—without fabricating a lid height. The series converges rapidly (Beychok notes that $m = 1, 2, 3$ are typically sufficient), and for the low effective stack heights in this system (9 m + plume rise, typically 15–25 m), the $g_3$ contribution is modest compared to $g_1 + g_2$ except when $L$ is very low (strong surface inversions), which is precisely when $S_{\text{inv}}$ is high and the empirical boost is largest. The empirical boost has been calibrated against field observations (Section 7.2) to produce correct odor-event predictions under inversion conditions.

C.2.4 Vertical Dispersion ($\sigma_z$) Not Distance-Dependent

Departure: The standard Pasquill-Gifford formulation computes $\sigma_z$ as a function of downwind distance $x$ and stability class, using power-law relationships $\sigma_z = c\,x^d + f$ with coefficients tabulated by Martin (1976) for classes A through F. At $x = 500$ m, $\sigma_z$ ranges from approximately 8 m (class F, very stable) to over 200 m (class A, very unstable). This system uses a single stability-dependent value $\sigma_z = 8 + 17(1 - s)$ evaluated at a representative distance of ~500 m, yielding 8 m (stable) to 25 m (unstable).

Criticism: The Pasquill-Gifford $\sigma_z$ grows continuously with distance, and this growth is the primary mechanism by which concentration decreases with downwind distance in the standard model (via $C \propto 1/\sigma_z$). By using a fixed $\sigma_z$, the model decouples vertical mixing from transport distance, losing the physical connection between plume age and vertical spread. The unstable value of 25 m is far below the Pasquill-Gifford class A value of ~223 m at 500 m, suggesting the model substantially underestimates vertical mixing in strongly convective conditions.

Justification: The fixed $\sigma_z$ is used only within the ground-level fraction $F_g$ (Section 5.4), not in the distance decay calculation. The distance-dependent dilution that $\sigma_z(x)$ provides in the standard model is handled separately by the BVF-adjusted decay exponent $\beta(s)$ (Section 4.3.3), which produces the correct qualitative behavior: fast decay under unstable conditions (both $\sigma_y$ and $\sigma_z$ growing, $\beta = 2.2$) and slow decay under stable conditions ($\sigma_z$ capped by BVF, $\beta = 1.0$). The fixed $\sigma_z$ in $F_g$ answers a different question: “at a typical community exposure distance, what fraction of the elevated plume reaches ground level?” For this purpose, a single representative value is sufficient.

The 25 m unstable ceiling reflects an intentional conservatism. Pasquill-Gifford class A ($\sigma_z \approx 223$ m at 500 m) represents extreme daytime convection with strong solar heating and light winds—conditions that are uncommon and that simultaneously produce strong buoyant plume rise (lifting emissions away from ground level) and vigorous mechanical dilution. Under such conditions, the system’s reduced $\sigma_z$ slightly overpredicts the ground-level fraction, erring on the side of caution. For a community health information tool, overprediction (false alerts) is preferable to underprediction (missed odor events).

C.2.5 Distance Decay: Power Law vs. $1/(\sigma_y \cdot \sigma_z)$

Departure: In the standard model, ground-level concentration decreases with downwind distance through the growth of $\sigma_y(x)$ and $\sigma_z(x)$, giving $C \propto 1/(\sigma_y \sigma_z)$. With Pasquill-Gifford power-law growth, this yields approximately $C \propto x^{-1.5}$ to $x^{-2}$ depending on stability. This system uses a truncated power law: $f_{\text{distance}} = (1 - d/d_{\max})^{\beta(s)}$ where $\beta(s) = 2.2 - 1.2\,s$ (Section 4.3.3). This function reaches zero at $d_{\max}$, whereas the standard Gaussian has an infinite tail.

Criticism: The truncated power law has a fundamentally different functional form from $1/(\sigma_y \sigma_z)$. The hard cutoff at $d_{\max}$ has no physical basis—Gaussian plumes have infinite (if vanishingly small) concentration at all distances. The decay curvature near the source and at long range differs from the standard model, potentially producing incorrect spatial gradients in the odor score.

Justification: The infinite Gaussian tail is physically unrealistic at distances beyond a few kilometers for the volatile organic compounds modeled here. At those distances, chemical transformation (photolysis, hydroxyl radical reaction), dry deposition, and gravitational settling—none of which the standard steady-state Gaussian model accounts for—reduce concentrations below detection. The $d_{\max}$ cutoff (5–10 km, stability-dependent) acts as a pragmatic proxy for these removal mechanisms. Within the truncated range, the exponent $\beta(s)$ correctly reproduces the key physics: stable conditions produce slower decay (BVF-capped $\sigma_z$, concentration falls as $\sim 1/x$) while unstable conditions produce faster decay (both $\sigma_y$ and $\sigma_z$ growing, concentration falls as $\sim 1/x^2$). The transition between these regimes is continuous and monotonic.

C.2.6 Lateral Dispersion: Angular $\sigma$ vs. $\sigma_y(x)$

Departure: The standard model uses $\sigma_y$ in meters, growing with downwind distance. This system uses a wind-speed-dependent angular spread parameter $\sigma$ in degrees (50° at $\leq$8 mph, 30° at 8–15 mph, 20° at $>$15 mph) applied as $f_{\text{angular}} = \exp[-0.5(\Delta\theta / \sigma)^2]$ (Section 4.3.2).

Criticism: The angular approach does not grow $\sigma_y$ with distance. In the standard model, a plume at 5 km downwind has a much wider crosswind spread (in meters) than at 500 m, even if the angular width is similar. The wind-speed dependence of $\sigma$ conflates two distinct effects: turbulent mixing (which determines $\sigma_y$) and mean-flow meandering (which is handled by the separate meander fraction). The three discrete wind-speed bins introduce discontinuities.

Justification: The angular formulation is equivalent to the standard $\sigma_y$ approach in the far field. For a receptor at distance $x$ and crosswind offset $y$, the angular offset is $\Delta\theta \approx y/x$ (small-angle approximation). Substituting into the angular Gaussian gives $\exp[-0.5(y/(x\sigma_\theta))^2]$, which is identical to the standard form $\exp[-y^2/(2\sigma_y^2)]$ with $\sigma_y = x \cdot \sigma_\theta$. That is, $\sigma_y$ does grow linearly with distance implicitly through the angular parameterization. The wind-speed dependence captures the well-documented relationship between wind speed and lateral plume spread: stronger winds produce narrower, more directional plumes. The 50°/30°/20° bins are empirically calibrated and consistent with Pasquill-Gifford angular spread measurements, which show lateral spread angles of 25°–60° for class A/B and 8°–15° for class E/F. The discrete bins are smoothed by the continuous Gaussian falloff at bin boundaries.

C.2.7 Briggs Plume Rise: Coefficient 2.66 vs. 2.6

Departure: The stable plume rise formula (Section 4.3.1) uses a coefficient of 2.66: $\Delta H = 2.66\,(F_b / (U \cdot N^2))^{1/3}$. Briggs (1969) originally published the coefficient as 2.6. Later refinements by Briggs (1975) and subsequent AERMOD documentation use values between 2.6 and 2.7.

Criticism: Using 2.66 instead of 2.6 overpredicts stable plume rise by approximately 2.3%, which increases the effective stack height and reduces predicted ground-level concentration. The choice of 2.66 is not standard and should be cited.

Justification: The difference between 2.6 and 2.66 produces a 2.3% change in $\Delta H$, which propagates through $\exp[-H_e^2/(2\sigma_z^2)]$ to an even smaller change in ground-level fraction (typically $< 1\%$ for the stack heights and conditions in this system). This is well within the uncertainty bounds of both the buoyancy flux estimation and the Brunt-Väisälä frequency parameterization. The coefficient 2.66 is taken from AERMOD’s implementation (EPA, 2024; prise.f source code, v24142), which represents the current regulatory standard. The variation across published sources (2.6–2.7) reflects the empirical nature of the original derivation.

C.2.8 Buoyancy Flux: $T_s$ vs. $T_a$ in Denominator

Departure: The buoyancy flux calculation (Section 4.3.1) uses: $F_b = g \cdot v_s \cdot (d/2)^2 \cdot (T_s - T_a) / T_s$ Some references (e.g., Briggs 1969, Turner 1994) define the buoyancy flux with $T_a$ (ambient temperature) in the denominator rather than $T_s$ (stack temperature).

Criticism: Using $T_s$ in the denominator instead of $T_a$ systematically underestimates the buoyancy flux. For a stack at 200°C ($T_s = 473$ K) and ambient 20°C ($T_a = 293$ K), the ratio is $F_b(T_s) / F_b(T_a) = T_a / T_s = 293/473 \approx 0.62$, a 38% reduction. This significantly reduces computed plume rise and therefore increases predicted ground-level concentration.

Justification: The AERMOD source code (prise.f, v24142) uses $T_s$ in the denominator, consistent with this implementation. The physical rationale is that the buoyancy of a parcel depends on the density deficit relative to ambient air, which for an ideal gas is $(T_s - T_a)/T_s$ when evaluated in the stack frame and $(T_s - T_a)/T_a$ when evaluated in the ambient frame. Both formulations appear in the literature, and the difference is partially compensated by the empirical coefficients (2.6–2.7) that were fit to observations using a specific form. Using $T_s$ with the AERMOD coefficient 2.66 is internally consistent. Moreover, the lower plume rise from the $T_s$ formulation produces higher predicted ground-level concentrations, which is conservative for a community health tool.

C.2.9 Log-Law Wind Profile: Missing Zero-Plane Displacement and Stability Correction

Departure: The log wind profile (Monin & Obukhov, 1954; Stull, 1988, ch. 9) in its full form is: $u_z = (u_*/\kappa)[\ln((z - d)/z_0) + \psi(z, z_0, L)]$ where $d$ is the zero-plane displacement (m), $\psi$ is the Monin-Obukhov stability correction, and $L$ is the Obukhov length. This system uses the neutral, zero-displacement form: $U(z) = U_{\text{ref}} \cdot \ln(z/z_0) / \ln(z_{\text{ref}}/z_0)$ with $z_0 = 0.5$ m and $z_{\text{ref}} = 10$ m (Section 4.3.2).

Criticism: Omitting the zero-plane displacement $d$ (typically 2/3 to 3/4 of obstacle height; Oke (1987) gives $d \approx 20$ m for a 30 m forest canopy) overestimates wind speed at low heights in areas with significant surface obstacles. For the suburban/light-industrial terrain around the Pineville facility, $d$ could be 3–5 m, which would reduce the wind speed ratio and therefore reduce estimated dilution.

Omitting the Monin-Obukhov stability correction $\psi$ is more significant. Under stable conditions (which are precisely the conditions most favorable for odor events), $\psi > 0$, which reduces the wind speed at plume height relative to the neutral assumption. The model therefore overestimates dilution during stable nighttime conditions—the opposite of conservative.

Justification: The effective stack height in this system is low: 9 m physical height plus typically 5–15 m of Briggs plume rise, yielding $H_e \approx 14$–24 m. At these heights, the zero-plane displacement effect is modest. Stull (1988) notes that the log-law is “generally considered to be a more reliable estimator of mean wind speed than the wind profile power law in the lowest 10–20 m,” which is precisely the relevant range. The extrapolation cap of 2.0 at line 444 of prediction.js prevents unrealistic results.

The missing stability correction is a genuine limitation. However, the effect is partially compensated by two other mechanisms: (1) the stability-dependent dilution efficiency $\eta(s) = 1.0 - 0.35\,s$ (Section 4.3.3), which reduces dilution under stable conditions, and (2) the inversion factor $I$ (Section 4.4), which directly amplifies scores during the stable nighttime conditions where the $\psi$ correction would be most important. These are empirical compensations rather than first-principles corrections, but they act in the correct direction and have been calibrated against field observations.

C.2.10 Surface Roughness Length ($z_0 = 0.5$ m)

Departure: The system uses a single roughness length $z_0 = 0.5$ m for the entire study area. Wieringa (1992) gives typical ranges of 0.1–0.5 m for suburban terrain and 0.5–1.0 m for brush/forest.

Criticism: A single $z_0$ cannot represent the heterogeneous terrain of the study area, which includes open industrial yards, forested bluffs, residential neighborhoods, and the Red River floodplain. Using $z_0 = 0.5$ m (upper end of suburban) may underestimate wind speed over open terrain (floodplain, $z_0 \approx 0.01$–0.05 m) and overestimate it over forested areas ($z_0 \approx 0.5$–1.0 m).

Justification: The KAEX ASOS station that provides reference wind measurements is located at an airfield with open fetch. The log-law ratio corrects from this 10 m measurement height to the plume effective height. Since the plume travels over mixed terrain between the source and receptor, the relevant $z_0$ is an effective average over the transport path. The value 0.5 m is consistent with EPA guidance for “urban” or “suburban” land use in AERMOD applications and represents a reasonable central estimate for the mixed Pineville–Alexandria landscape. A spatially varying $z_0$ field would require a land-use database and sector-dependent roughness calculation, adding complexity disproportionate to the improvement for a relative scoring model.

C.2.11 AERMOD Meander: Default Timescale and $\sigma_v$ Parameterization

Departure: The meander fraction calculation (Section 4.3.2) reproduces the AERMOD MEANDR subroutine with a meander timescale $T = 24$ hours and lateral turbulence $\sigma_v = 0.5 \cdot \min(U,\, 2.0)$ m/s. The EPA’s LOWWIND white papers (2021) identified issues with the default AERMOD meander formulation and introduced alternatives: LOWWIND1 uses a fixed minimum $\sigma_v = 0.5$ m/s, LOWWIND2 reduces the timescale to 12 hours and uses minimum $\sigma_v = 0.3$ m/s, and both cap the FRAN upper limit at 0.95 instead of 1.0.

Criticism: The default $\sigma_v = 0.2$ m/s minimum in AERMOD was found to be “too low by at least a factor of 2” (EPA, 2021). This system’s $\sigma_v$ formula yields a minimum of ~0.25 m/s (at $U = 0.5$ m/s floor), which is closer to the LOWWIND2 minimum but still below the Hanna (1990) observational finding that $\sigma_v$ maintains approximately 0.5 m/s even as wind speed approaches zero. The 24-hour timescale may overestimate meander at intermediate distances compared to the LOWWIND2 12-hour value. The FRAN upper limit is not capped at 0.95.

Justification: This system’s $\sigma_v$ parameterization is wind-speed-dependent rather than using a fixed minimum, which better captures the physical relationship between mean flow and turbulent fluctuations. At the low wind speeds where meander is most important ($U < 2$ m/s), the formula produces $\sigma_v = 0.25$–0.5 m/s, which falls between the AERMOD default and LOWWIND1 values. The 24-hour timescale matches the AERMOD default and is appropriate for the spatial scales of this application (receptors at 0.5–5 km). The FRAN upper limit of 1.0 vs. 0.95 produces a negligible difference in practice, as it only affects the transition to fully circular dispersion under near-calm conditions where the score is already dominated by the circular pooling term. Adopting the LOWWIND2 timescale of 12 hours would modestly improve low-wind performance and could be implemented as a future refinement.

C.2.12 Brunt-Väisälä Frequency: Parameterized vs. Measured

Departure: AERMOD computes the Brunt-Väisälä frequency $N$ from the vertical potential temperature gradient in the upper-air sounding: $N = \sqrt{(g/\theta)(\partial\theta/\partial z)}$. This system parameterizes $N$ as a linear function of the stability parameter: $N = \max(0.005,\, s \times 0.04)$ s$^{-1}$, where $s \in [0, 1]$ is itself derived from surface observations (Section 4.3.3).

Criticism: The BVF is a property of the vertical temperature profile, not the surface layer. Surface-based proxies (temperature trend, dewpoint depression, wind speed) cannot fully characterize the stratification aloft. The linear mapping from $s$ to $N$ is a crude approximation: real BVF values can vary from near-zero (well-mixed boundary layer) to $>$0.04 s$^{-1}$ (strong radiation inversions), with complex vertical structure (e.g., elevated inversions with neutral layers below).

Justification: Direct BVF measurement requires upper-air soundings (radiosondes), which are not available in real time at this location. The nearest radiosonde site (Shreveport, KSHV) is ~200 km away and launches only twice daily (00Z and 12Z), providing inadequate spatial and temporal resolution for hourly community-level predictions. The parameterized $N$ captures the first-order relationship between surface stability indicators and vertical stratification. The range 0.005–0.04 s$^{-1}$ spans the physically plausible range for the lower boundary layer in central Louisiana’s subtropical climate. The linear mapping is conservative: it avoids the extreme values that could produce unrealistically long or short plume transport distances. Full AERMOD implementations also face BVF estimation challenges when sounding data is unavailable, and EPA guidance permits the use of surface-based stability estimates (AERMET preprocessor) as surrogates.

C.3 Equations Verified Against Published Sources

The following equations were verified term-by-term against published reference material and found to match exactly.

EquationThis SystemReferenceStatus
Buoyancy flux $F_b = g\, v_s (d/2)^2 (T_s - T_a)/T_s$ AERMOD prise.f v24142 Exact match
Stable plume rise $\Delta H = 2.66\,(F_b / (U \cdot N^2))^{1/3}$ Briggs (1975); AERMOD implementation Exact match
Unstable plume rise $\Delta H = 1.6\, F_b^{1/3}\, x_f^{2/3} / U$ Briggs (1969, 1975) Exact match
Distance to final rise ($F_b < 55$) $x_f = 49\, F_b^{5/8}$ Briggs (1969): $3.5 \times 14 \times F_b^{5/8}$ Exact match
Distance to final rise ($F_b \geq 55$) $x_f = 119\, F_b^{2/5}$ Briggs (1969): $3.5 \times 34 \times F_b^{2/5}$ Exact match
Log-law wind ratio $U(z_2) = U(z_1) \cdot \ln(z_2/z_0) / \ln(z_1/z_0)$ Stull (1988), eq. 9.7.5; Monin & Obukhov (1954) Exact match (with $d = 0$)
Ground-level Gaussian $F_g = \exp[-H_e^2/(2\sigma_z^2)]$ $g_1$ at $z = 0$ (Pasquill, 1961; Turner, 1994) Exact match (see C.2.2)
Crosswind Gaussian $f = \exp[-0.5(\Delta\theta/\sigma)^2]$ Equivalent to $\exp[-y^2/(2\sigma_y^2)]$ (see C.2.6) Equivalent form
AERMOD FRAN (meander) $(2\sigma_v^2 + U_{\text{mean}}^2(1 - e^{-t/T})) / U^2$ EPA AERMOD Low Wind White Paper (2021) Exact match

C.4 Summary of Conservatism Profile

A critical question for any community health information tool is whether systematic errors produce false negatives (missed odor events) or false positives (false alerts). The following table summarizes the direction of bias introduced by each departure.

DepartureEffect on Predicted OdorConservative?
Missing $g_2$ reflection factor Underestimates by constant factor Neutral (absorbed in calibration)
Empirical $g_3$ inversion boost Direction depends on conditions Calibrated to field data
$\sigma_z$ capped at 25 m (unstable) Overpredicts ground-level fraction Conservative (more odor predicted)
Truncated distance decay at $d_{\max}$ Underpredicts at extreme range Anti-conservative at $>d_{\max}$, but concentrations there are sub-threshold
$T_s$ in buoyancy flux denominator Less plume rise → more ground-level odor Conservative
No stability correction $\psi$ in log-law Overestimates dilution in stable conditions Anti-conservative (partially offset by $\eta(s)$ and inversion factor)
No zero-plane displacement $d$ Slightly overestimates wind at low heights Slightly anti-conservative
Parameterized BVF instead of measured Direction varies; smooths extremes Neutral (avoids both over- and under-prediction of rare extremes)

The conservatism profile is mixed but leans toward overprediction. The two anti-conservative departures (missing log-law stability correction and zero-plane displacement) are both in the wind profile, which affects dilution. These are partially compensated by the stability-dependent dilution efficiency and inversion amplification factors. The remaining departures are either conservative (reduced $\sigma_z$, $T_s$ in buoyancy flux), neutral (absorbed in calibration), or negligible in magnitude (Briggs coefficient, FRAN cap). On balance, the model is more likely to alert when no odor occurs than to miss an odor event—an appropriate bias for a community information tool.

C.5 Literature Cited in This Appendix

Version History

Version Date Changes
1.0 2026-03-09 Initial release. Basic Gaussian plume model with temperature-driven emission factors, single-source wind transport, and static time-of-day modifiers.
1.1 2026-03-13 Added radiation inversion detection using Brutsaert radiative cooling physics and METAR cloud ceiling extraction. Introduced inversion type classification (radiation, fog, subsidence, residual, surface). First field calibration against observed odor events.
1.2 2026-03-15 Rebalanced Pineville nighttime emission model using LDEQ production data (Incident #227489, 2026-01-28). Ground-level cylinder operations elevated from secondary to co-dominant nighttime source.
1.3 2026-03-19 Added Open-Meteo ECMWF forecast blending to supplement NWS data. Time-dependent blend weights favor ECMWF at longer lead times. Added community odor observation logger for model calibration. Solar position model replaced hardcoded daytime hours with Spencer (1971) astronomical calculations.
1.4 2026-03-23 Data quality fixes for observation pipeline. Calm wind (0 mph) no longer stored as “north” — direction nulled so dispersion model uses circular pooling mode. Odor submissions now fetch live KAEX METAR instead of relying on 30-min-stale CRON data. Removed strftime wrappers from D1 range queries to enable index usage.