An adaptation of the monte carlo plume model described by Alfred K Blackadar.
Introduction
The DIRSIG Blackadar plume model is an enhanced implementation of
the MonteCarlo smoke plume model described in the book "Turbulence
and Diffusion in the Atmosphere" by Alfred K. Blackadar. In Appendix
E, Blackadar describes "A Monte Carlo Smoke Plume Simulation" and
the book includes a simple program written in
BASIC that can model a smoke
plume as a series of particles that are released and perform a
"drunkards walk". All equations and some descriptive text included
in this document have been derived or extracted from Blackadar’s
book and are noted with references. Blackadar’s published BASIC
code is the basis for the enhanced C++ implementation manifested
in the blackadar
utility included in DIRSIG. The primary addition
in our enchanced model to Blackadar’s particle (zero volume) approach
is that we are modeling puffs (nonzero volume).
History
The original implementation of the Blackadar plume model was included
in DIRSIG4. In the DIRSIG4 era, the plume was internally ngenerated
at runtime and internally modeled as a series of puffs. In the DIRSIG5
era, the plume is externally generated (by the blackadar
command
line program), stored in an OpenVDB file and
then read in by a plume specific plugin.
Theory
Each puff is a sphere that starts out with a specific diameter. The motion of each puff uses the same particle motion that was described by Blackadar. However, as the puff travels dispersion (from turbulence in the atmosphere) causes the puff to expand in volume. The rate of expansion is driven by dispersion coefficients which are usually different in the lateral and vertical directions. For the time being, the lateral and vertical dispersion coefficients are equal and the vertical dispersion coefficient is computed as a function of downwind range and the atmospheric stability.
When we model a turbulent wind field, there are two perspectives to consider. The first is the Eulerian perspective where we look at the statistics of the wind flow at a specific location. In the presence of a turbulent wind field, we expect the wind velocity to vary at a specific location. The second is the Lagrangian perspective where we look at the statistics of the wind field acting upon a point (particle) that moves with the flow. In the presence of a turbulent wind field, a particle traveling in that field will experience varying wind velocities as a function of time and, hence, location. In this model, we must consider an Eulerian modeling approach to describe the variation of the wind field at the stack and a Lagrangian modeling approach to describe the variation within the wind field as each particle/puff travels downwind.
As it was described earlier, the puff motion can be described as a "drunkards walk" which is described by traditional Lagrangian motion. A "random walk" would be where a particle "chooses" a totally random direction (uncorrelated to its previous direction) at each time step. In a drunkards walk, the particle "remembers" part of its previous direction (correlated) and adds a random (uncorrelated) component to the previous direction to determine the new direction. In this treatment, the puff velocity for the next time step is correlated to the current time step by combining a portion of the current velocity with a random velocity component.
\begin{equation} u( t + \Delta t) = W \cdot u(t) + u_{\mathrm{random}}(t) \end{equation}
where W is a weighting factor for the portion of the current velocity component. In the context of the "drunkards walk", this weighting factor is the fraction of the previous velocity that is remembered. The rate of memory loss is goverened by the Lagrangian time scale which is a function of the atmospheric stability classification; stable conditions will have the longest times scales and produce slow changes that result in the largest and longer lasting loops. Specifically, this weighting factor is defined by the Lagrangian autocorrelation function, \(R(\Delta t)\). The autocorrelation function describes the similarity as a function of time scales. For very small time scales we expect R → 1 (highly correlated) and for very long time scales we expect R → 0 (highly uncorrelated). For the X dimension, the velocity for the next time step can be written as:
\begin{equation} u_{x}( t + \Delta t) = u_{x}(t) R_{x}(\Delta t) + u_{x}'(t) \end{equation}
The random velocity component (\(u_x'(t)\)) is assumed to be Gaussian distributed with a zero mean and with a standard deviation (\(\sigma_{x}'\)) defined as:
\begin{equation} \sigma_{x}' = \sigma_{x} \sqrt{ 1  R_{x}^2( \Delta t) } \end{equation}
Note, that the variance for the random velocity is also correlated through the inclusion of the autocorrelation function. The Lagrangian autocorrelation function, \(R(\Delta t)\) is assumed to have the following form:
\begin{equation} R_{x}( \Delta t) = \exp \left ( \frac{\Delta t}{T_L} \right ) \end{equation}
where \(T_L\) is the Lagrangian integral time scale of the velocity. In this model, the Lagrangian integral time scale is computed directly from the stability of the atmosphere. For stable atmospheres, the time scale is long and for unstable atmospheres the time scale is short.
Model Usage
The blackadar
utility program is used to generate voxelized plumes
suitable for ingest by the PlumeVDB plugin.
The plume is described by a small XML document contraining a series
of variables. In DIRSIG4, this was stored inside the .scene
file.
In DIRSIG5, this can be stored in separate file that is provided
to the standalone blackadar
program.
<blackadarplume> <directinsolation> 600.0</directinsolation> <airtemperature> 300.0</airtemperature> <windspeed> 4.0</windspeed> <winddirection> 90.0</winddirection> <surfaceroughness> 0.01</surfaceroughness> <mixedaltitude> 1000.0</mixedaltitude> <stacklocation><point><x>0.0</x><y>0.0</y><z>0.0</z></point></stacklocation> <stackheight> 80.0</stackheight> <stackdiameter> 2.0</stackdiameter> <releasebuoyancy> 10.0</releasebuoyancy> <releasetemperature> 350.0</releasetemperature> <releaserate> 100.0</releaserate> <puffsperstep> 8</puffsperstep> <stepcount> 320</stepcount> <debugfilename>blackadar.txt</debugfilename> </blackadarplume>
The table following describes these variables.
Input Variable  Description  Suggested Value(s) 


The broadband, incident radiative flux [W/m^{2}] 
300.0  600.0 

The ambient air temperature [K] 
280.0  320.0 

The wind speed at the base of the stack [m/s] 
1.0  5.0 

The wind direction at the base of the stack [degrees, CW from N] 
0.0  360.0 

The surface roughness length [m] 
0.01  0.10 

The altitude of the mixed layer [m] 
500.0  2000.0 

The Scene ENU X/Y/Z location of the stack base [m] 
(any valid location) 

The height of the stack [m] 
60.0  150.0 

The diameter of the stack [m] 
1.0  4.0 

The Briggs buoyancy flux parameter (F) at the stack [m^{4}/s^{3}] 
4.0  8.0 

The exiting stack temperature [K] 
300.0  400.0 

The release rate of the exiting material [ppm/s] 
100.0  10,000 

The number of new puffs released each time step 
(user defined) 

The number of time steps to establish the plume 
(user defined) 
Notes

With DIRSIG5, the plugin instantiates the plume into the scene so it is recommented to leave the
stack_location
at the Scene ENU origin (0
,0
,0
) and set the origin in the PlumeVDB plugin. 
Since the plume is currently advanced in 1 second intervals. the
stepcount
can be interpreted as how many seconds have elapsed since the plume started to be emitted at the stack. 
The concentration of given puff is initially computed from the
releaserate
,puffsperstep
andsteptime
. The release rate is a concentration rate [ppm/second]. If we multiply thereleaserate
with thesteptime
(default =1
second) we would arrive at the concentration released per time step. Since multiple puffs are released in each time step, we must divide by thepuffsperstep
to obtain the concentration per puff. 
A higher value of
puffsperstep
will make the plume appear to be more geometrically continuous. For example, in a high wind condition, the user might want to increasepuffsperstep
so that the plume does not feature an artifical break. Think of it from a sampling perspective where thepuffsperstep
controls how many elements are used to represent a nearly continuous plume. If this number is too small, a break in the plume might simply be a result of not enough puffs were available to span the a distance. 
There is only one material in the plume at this time.

The Briggs buoyancy flux parameter (\(F_b\)) can be computed in various ways. The most practical is from the following equation:
\begin{equation} F_b = g \frac{V_s}{\pi} \frac{T_s  T_a}{T_s} = g v_s r^2 \frac{T_s  T_a}{T_s} \end{equation} where \(g\) is the gravity constant, \(V_s\) is the stack volume gas flow rate [m^3/s], \(v_s\) is the stack gas exit velocity [m/s], \(r\) is the stack radius [m], \(T_s\) is the stack gas temperature [K] and \(T_a\) is the ambient air temperature [K].
The blackadar
utility is run by supplying the name of the input
XML file and the name of the output VDB file:
blackadar
utility.$ blackadar plume.xml plume.vdb XML [plume.xml] parsed without errors Writing plume to: plume.vdb ... done
Since this plume model has a random component, there is an optional commandline argument to specify the random seed:
blackadar
utility with a random seed.$ blackadar random_seed=1 plume.xml plume.vdb Using random seed = 1 XML [plume.xml] parsed without errors Writing plume to: plume.vdb ... done
Model Assumptions
There are several assumptions and caveats that should be considered when using this plume model.

Presently, the plume is static during the simulation. The user supplies the number of time steps that will be used to evolve the plume during startup and then the plume is frozen during the remainder of the simulation. This is less than ideal because systems that build the image over time should see the plume changing with time. An improvement to this will be addressed at a future time.

The plume is unaware of the rest of the scene. As a result, the plume will not flow around other scene elements (e.g. buildings, terrain, etc.). Instead, it will evolve as if the world were perfectly flat and open. Hence, the plume will blow through a building as if it wasn’t there. As long as nothing obscures the view of the plume, the sensor will see it. In the case where a plume passes into a building, the building will block the view of the plume.

The release temperature is used purely for the radiometry side of the model. Changing the plume release temperature will not affect the buoyancy of the gas (as it should). At a future time, this release temperature and the gas buoyancy will be coupled.

The dispersion coefficient is currently assumed to be the same in the lateral and vertical dimensions. However, it is well known that the dispersion is different in the lateral and vertical dimensions.

At this time, you can only have one gas in the plume.

At this time, DIRSIG4 can only have one plume in the scene. In DIRSIG5, you can instantiate multiple instances of the plume but they do not interact or mix with each other.
Future Improvements
The following is a list of future improvements to be made to the model.

Convert some of the variables that are currently usersupplied to be automatically pulled from the weather descriptions that DIRSIG already has been supplied:

The
insolation

The
windspeed

The
airtemperature


Make the delta time (see the
steptime
variable) userdefined. 
Make the plume dynamic with time so that it changes during a sequence of captures (currently the plume evolves during initialization and is frozen during image generation).

Provide a way to specify the atmospheric stability class rather than the surface roughness (see the
surfaceroughness
variable), since the later is less intuitive to most users. 
Compute the Brigg’s buoyancy flux parameter (\(F_b\)) directly from the other parameters rather than making the user supply it. We would need to add either the stack volume flow rate or exit velocity as an input, but either of these would be more intuitive than the Briggs parameter.

Allow the lateral and vertical dispersion to be independent. The challenge there will be figuring out how to model a puff that is no longer a sphere.

Allow for multiple gases in a single plume.