An adaptation of the monte carlo plume model described by Alfred K Blackadar.


The DIRSIG Blackadar plume model is an enhanced implementation of the Monte-Carlo 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 (non-zero volume).


The original implementation of the Blackadar plume model was included in DIRSIG4. In the DIRSIG4 era, the plume was internally ngenerated at run-time 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.


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.

puff plot
Figure 1. A plot of an example plume showing the puffs as a series of expanding spheres.

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 stand-alone blackadar program.

Example XML document with the model input variables.
  <directinsolation>   600.0</directinsolation>
  <airtemperature>     300.0</airtemperature>
  <windspeed>            4.0</windspeed>
  <winddirection>       90.0</winddirection>
  <surfaceroughness>    0.01</surfaceroughness>
  <mixedaltitude>     1000.0</mixedaltitude>
  <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>

The table following describes these variables.

Table 1. Blackadar plume model parameters.
Input Variable Description Suggested Value(s)


The broad-band, incident radiative flux [W/m2]

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 [m4/s3]

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)


  • 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 and steptime. The release rate is a concentration rate [ppm/second]. If we multiply the releaserate with the steptime (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 the puffsperstep 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 increase puffsperstep so that the plume does not feature an artifical break. Think of it from a sampling perspective where the puffsperstep 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:

Example usage of the 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 command-line argument to specify the random seed:

Example usage of the 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 user-supplied 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) user-defined.

  • 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.