Prepared 5 September 2026. Engine version 0.2.0.
The second prototype uses measured-source elevation data to construct a three-dimensional air grid, redirects an assumed wind field around and over that terrain, and transports separate gas and particle tracers through the resulting air cells. It tracks emitted, airborne, deposited and escaped mass. It can retain a plume after a flare stops and integrate concentrations over full one-hour and 24-hour periods.
This is an exploratory numerical model, not a calibrated Woodfibre exposure prediction. Its wind calculation conserves air volume but does not solve atmospheric momentum, temperature, pressure systems or turbulence. Its inversion is an imposed restriction on vertical exchange, not a predicted temperature inversion. Near-surface output is the mean within the lowest model air cell, not a measurement or prediction at breathing height. Coarse-grid numerical spreading is a material limitation.
The terrain input is Natural Resources Canada's legacy CDEM tile 092G. The archived raster has 0.75-arcsecond horizontal spacing, NAD83(CSRS) horizontal coordinates and CGVD28 elevations. The prepared 20 × 20 km extract contains 100 × 100 cell-centred, 200 m representative elevations in a local azimuthal-equidistant projection centred at 49.703° N, 123.22° W. Its elevations range from 0 to 2,166.1 m. Each representative elevation averages 25 bilinear samples from the native raster; optional further model coarsening uses area-weighted means. This preserves the source's broad terrain but smooths shoreline cliffs, small valleys and structures. CDEM is historical and is not a survey of the completed facility. NRCan data catalogue, official tile download.
For each horizontal column, the model intersects fixed elevation layers with that column's terrain height. An air cell extends from
z_low = max(terrain elevation, nominal layer bottom)
z_high = nominal layer top
volume = dx * dy * (z_high - z_low)
Intervals entirely inside the resolved ground are absent. A bottom interval thinner than 8 m is merged with the cell above it to avoid arbitrarily small time steps. Ground heights are not lowered to create passageways. The resulting bottom air-cell depth varies by location and must be displayed with concentration values. On high slopes this depth can be hundreds of metres; calling those values “2 m concentrations” would be incorrect.
Neighbouring columns exchange air and pollutant only through the vertical portions of their cell faces that overlap above both terrain surfaces. Rock portions have no faces and therefore carry no flow. Vertical connections occur between adjacent cells in each column. The model top lies above every mountain; the default altitude is at least 3,600 m. Surface topography is represented as vertical steps between averaged columns, rather than exact sloping faces. The fine source raster is not the simulation's effective resolution.
The user supplies a wind direction from which the air comes, measured clockwise from north, and a speed. The initial horizontal vector is:
u0 = -speed * sin(direction_from)
v0 = -speed * cos(direction_from)
w0 = 0
Separate upper-level direction and speed can be supplied above a selected inversion altitude. The optional transition depth blends the vectors. These are prescribed scenario conditions, not local observations reproduced throughout the domain.
For each internal open face, the initial volume flux is F0 = normal wind × face area, in m³/s. The engine seeks a correction that balances the sum of inflows and outflows in every air cell. With a pressure-like mathematical potential p, the corrected flux from cell i to cell j is:
Fij = F0ij + gij * (pi - pj)
gij = mobility * face_area / centre_distance
sum(outward corrected fluxes for cell i) = 0
There is no flow through terrain. Lateral boundary normal velocities retain their prescribed values. The top has zero potential, permitting compensating vertical volume flux. A conjugate-gradient solve with a vertical-line preconditioner produces the potential. The engine reports residual, volume-balance error and maximum face speed.
Changing the DEM changes the open air passages and pressure correction, so terrain actually affects the computed flow. This is not a Gaussian plume bent along a drawn valley. However, satisfying continuity alone does not establish the correct wind speed, direction, sheltering or circulation in Howe Sound. For example, an examined 3/5 m/s boundary-wind case produced local face speeds around 21 m/s. Such amplification is an uncalibrated consequence of this idealized boundary problem, not evidence that those local winds occur.
Mass-consistent diagnostic downscaling is an established model category; the Forest Service's WindNinja describes prescribing an initial field and adjusting it over terrain. This prototype does not run WindNinja or inherit its validation, surface treatment or thermal parameterizations. A recent primary evaluation also found that mass-conserving downscaling can reduce skill in thermally driven downslope events, underlining the difference between mechanical terrain adjustment and full atmospheric dynamics. WindNinja method overview, Seto and colleagues, 2025 evaluation.
An inversion altitude is specified in metres above sea level, not metres above each hillside. The engine reduces vertical correction mobility at and above that altitude by inversionStrength; a value of 0.05 means 5% of the unrestricted vertical mobility remains. The separate inversionMixingFraction controls remaining vertical pollutant diffusivity; 0.03 means 3% remains. These dimensionless values are sensitivity parameters without site calibration.
The user can prescribe a lower outflow below the cap and a different upper wind above it, allowing an elevated injected plume to encounter different transport from near-surface air. The default cap does not imply that a real inversion exists at that height. Inversions intersecting thin merged bottom cells are resolved only to the grid's vertical discretization.
The model does not calculate surface cooling, drainage generation, sea-breeze development, internal gravity waves, inversion breakup, hydraulic transitions, mountain-wave acceleration or lee-side turbulence. Outflow and inflow are forced boundary scenarios. Under zero wind the projected flow remains zero; diffusion and particle deposition can still operate. This allows accumulation in the numerical model, but it does not reproduce real weak-wind meandering or local thermal circulations.
Howe Sound field research supports resolving vertical and along-fjord variation: Jackson and Steyn's observations documented substantial spatial and vertical changes in outflow. That evidence motivates layers; it does not supply universal wind profiles or inversion settings for this prototype. Jackson and Steyn, 1994, Howe Sound field observations.
The state variable is pollutant mass Mi in grams within each air cell. Cell concentration is Ci = Mi / Vi, in g/m³. Each internal face transfers:
transfer_ij = dt * [Fij * C_upwind + Dij * (Ci - Cj)]
Dij = diffusivity * face_area / centre_distance
The identical mass is subtracted from one cell and added to the other. Horizontal and vertical diffusivities are supplied independently; example values of 10 and 1 m²/s are assumed, not inferred from weather or fitted to monitoring. Vertical diffusivity is multiplied by the inversion mixing fraction at and above the specified cap.
The model uses first-order upwind advection and explicit face diffusion. It selects a time step no larger than 0.8 times each cell's volume divided by the sum of its outgoing advective, diffusive and deposition coefficients. Incoming boundary air contains zero project tracer; outflow carries the cell's tracer concentration out of the domain. No external diffusion is imposed. A release at rate Q g/s adds exactly Q × active time to one specified air cell. The release stops at the chosen duration, while existing tracer keeps moving.
The mass identity is:
initial + emitted = airborne + escaped + deposited
The engine reports floating-point residual and any negative-mass clipping separately. Normal verified runs needed no clipping. The top and map sides are open; escaped mass is outside the model, not chemically destroyed. Domain-edge escape means later recirculation into the study area is not represented.
Numerical spreading matters: first-order upwind transport smears a sharp plume. Its effective numerical diffusion scales approximately with wind speed times grid width; at several hundred metres grid spacing it can exceed the chosen physical horizontal diffusivity. Numerical conservation does not eliminate that error. Concentration contour locations and community overlaps should therefore be treated as scenario illustrations. A grid-refinement and higher-order transport study is necessary before relying on detailed reach or peak concentrations.
The injected height is an assumption about the effective release location above the resolved source ground, or an explicitly supplied altitude. It is not calculated from flare heat. Injection is distributed over the selected air cell, so the actual lower and upper cell elevations are reported. Changing pollutant emission rate leaves wind, heat and injection geometry fixed. This provides valid linear scaling within the model but does not simulate how a larger real flare changes buoyant rise.
Gas channels are passive tracers without chemical conversion. CO, NOx, SO2 or a named VOC therefore share a unit-mass transport solution when assigned identical deposition parameters. This must not be described as chemical speciation or a chemical forecast.
The particle channel can add a prescribed downward settling speed on vertical faces and a surface deposition velocity. The lowest cell loses (settling speed + deposition velocity) × horizontal ground area × concentration to the deposited-mass ledger. Zero settling and deposition exactly recover the passive gas result. A value such as 0.001 m/s is a sensitivity assumption, not a measured effective PM2.5 settling speed. Particle size distributions, coagulation, humidity growth, rain scavenging and secondary aerosol formation are absent. Deposition on vertical terrain steps, vegetation and buildings is not separately resolved.
The selected elapsed-time map displays the state at that time. It is not an hourly mean. If averaging is requested, the model continues to the end of every full window, including the entire plume tail after the release stops. It does not divide a short simulation by 24 hours while assuming the missing hours are pollutant-free.
Concentrations are integrated using the trapezoidal rule over the resolved time steps. Time steps end exactly at release-stop, requested snapshot, hourly-block and averaging-window boundaries. Requested cumulative averages are:
mean_C(0,T) = integral(C(t), t=0..T) / T
When 24 hourly blocks are requested, the engine also reports the maximum of the means for 0–1 h, 1–2 h, …, 23–24 h, and the start of the maximum block for each cell. A map of these maxima combines different times at different locations; it is not one simultaneous plume. It is not a rolling one-hour maximum and cannot establish annual-percentile or multi-year standard attainment. All 24 hours use the same imposed scenario unless a future version supplies time-varying meteorology.
Objective overlays must use the corresponding pollutant and averaging period. A passive NOx mass field must not silently become NO2; a full-conversion screen requires an explicit NO2-equivalent emission mass basis and conversion label. Background air pollution is absent. A source increment below an objective does not establish that total ambient concentration is below it, and a model overlap does not establish an individual's dose or probability of harm. B.C. publishes specific guidance on NOx-to-NO2 conversion and the use of refined models for complex coastal flow. B.C. NO2 modelling guidance, B.C. dispersion guideline, July 2026.
Numerical checks cover:
An actual-DEM 48 × 48 test used 20,858 air cells and 58,474 internal faces. A 30-minute release of 1 g/s produced exactly 1,800 g. In the selected test's 24-hour integration, approximately 1,786.990 g escaped, 13.010 g deposited and negligible mass remained; relative mass-accounting error was approximately 4.4 × 10⁻¹². These numbers verify that particular numerical run, not Woodfibre's emissions or environmental effects. A tested inversion parameter change reduced mass entering the below-cap region from 169.8 g to 3.53 g; this verifies the implemented mechanism, not the strength of a local inversion.
A 48-to-64-column refinement check did not show convergence. For the same imposed capped 3/5 m/s wind case and 1 g/s passive source, the one-hour snapshot's maximum near-surface concentration increased from 0.543 to 0.838 µg/m³, while the Squamish marker decreased from approximately 0.00123 to 0.000671 µg/m³. Maximum diagnostic face wind increased from 21.2 to 27.2 m/s. Changes in represented terrain, cut-cell geometry, source placement and numerical spreading all contribute. These differences preclude treating the current contour positions or values as resolved exposure predictions.
The next scientific step is representative three-dimensional wind and turbulence data, verified source time histories and heat release, appropriate background/chemistry, larger meteorological coverage, spatial and time-step convergence, and comparison against independent wind and pollutant observations. A full model should reproduce observed inflow/outflow profiles and reversal timing before its concentration maps are used for public exposure conclusions. WRF/CALMET/CALPUFF or an appropriately evaluated momentum-resolving alternative remains the route to an assessed site model; this prototype provides inspectable mechanics and a data framework, not a substitute for that work.
The standalone engine is available as the browser/worker global WoodfibreFV or a CommonJS module. Objects use plain records and typed arrays and can be passed through worker messages.
const terrain = WoodfibreFV.resampleTerrain(dem, 48, 48);
const grid = WoodfibreFV.createGrid({
nx: 48, ny: 48, dx: 20000 / 48, dy: 20000 / 48,
elevations: terrain, inversionHeight: 500,
source: { x: -2467.88, y: -4197.82 }
});
const flow = WoodfibreFV.prepareFlow(grid, {
directionFrom: 200, speed: 3,
upperDirectionFrom: 225, upperSpeed: 5,
inversionHeight: 500, inversionStrength: 0.05
});
const result = WoodfibreFV.simulate(grid, flow, {
elapsedSeconds: 3600, releaseSeconds: 1800,
sourceHeightAGL: 200, emissionRateGs: 1,
averagingWindowsSeconds: [3600, 86400], hourlyBlockCount: 24,
species: { kind: 'pm', settlingVelocity: 0.001, depositionVelocity: 0 }
});
result.groundLayerMeanUgM3 and result.concentrationUgM3 describe the selected elapsed-time snapshot. result.averages holds full-window means. result.groundMaxHourlyUgM3 and result.groundMaxHourlyStartSeconds describe the hourly composite. result.integrationMass holds the full-horizon mass ledger; top-level mass fields describe the selected snapshot. grid.cells.zLo/zHi/volume and grid.columnNodes support vertical sections. sampleAt and section report cell means without pretending to interpolate to breathing height.