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 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 enhanced model to Blackadar’s particle (zero volume) approach is that we are modeling puffs (non-zero volume).

History

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

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 governed 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 tool in the DIRSIG vdb_tool utility program is used to generate voxelized plumes suitable for ingest by the PlumeVDB plugin. The plume is described by a small XML document containing 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.
<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 following table describes these variables.

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

insolation

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

300.0 - 600.0

airtemperature

The ambient air temperature [K]

280.0 - 320.0

windspeed

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

1.0 - 5.0

winddirection

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

0.0 - 360.0

surfaceroughness

The surface roughness length [m]

0.01 - 0.10

mixedaltitude

The altitude of the mixed layer [m]

500.0 - 2000.0

stacklocation

The Scene ENU X/Y/Z location of the stack base [m]

(any valid location)

stackheight

The height of the stack [m]

60.0 - 150.0

stackdiameter

The diameter of the stack [m]

1.0 - 4.0

releasebouyancy

The Briggs buoyancy flux parameter (F) at the stack [m4/s3]

4.0 - 8.0

releasetemperature

The exiting stack temperature [K]

300.0 - 400.0

releaserate

The release rate of the exiting material [ppm/s]

100.0 - 10,000

puffperstep

The number of new puffs released each time step

(user defined)

stepcount

The number of time steps to establish the plume

(user defined)

emitterfilename

The name of an ASCII/Text file containing emitter info

(user defined)

debugfilename

The name of an ASCII/Text file containing debug info

(user defined)

Discrete vs. Extended Source

The model was originally written to model a single source emitter location (e.g., a factory stack), where the puffs are emitted from a small region defined by the <stacklocation>, <stackheight> and <stackdiameter> (puffs are randomly generated within the stack opening).

plume point
Figure 2. A plume (in Blender) from a single emitter.

The model was extended in 2024 to support an arbitrary set of source emitter locations in order to model an extended source (i.e., the smoke plume from the fire line in a forest/wild fire). The extended source option is activated by setting the <emitterfilename> element to the name of an ASCII/Text file that contains the 3D locations of an arbitrary number of emitters. Each line in the file contains the 3D location (in the Scene ENU coordinate system) and an optional probability that it generates a puff at any given time. The puffs are generated using the following method:

  • The <puffsperstep> variable still drives how many puffs are released in each time step. The puffs generated in each time step are distributed to start from the emitters using a uniform distribution. Once a random emitter location has been selected, a second random number is draw to check against the probability for that emitter. If the emitter does not emit for the current time step, the process is repeated to select another candidate emitter.

  • The <stackdiameter> is still used to define the initial puff radius.

The example file below shows a set of 21 point emitters that are distributed along the Y-axis. The optional 4th token is not present for these emitters, so the probability is assumed to be 1 for each of them.

An example extended source emitter file containing 21 locations along a straight line and without (equal) probabilities.
0.0 -10.0 5.0
0.0  -9.0 5.0
0.0  -8.0 5.0
0.0  -7.0 5.0
0.0  -6.0 5.0
0.0  -5.0 5.0
0.0  -4.0 5.0
0.0  -3.0 5.0
0.0  -2.0 5.0
0.0  -1.0 5.0
0.0   0.0 5.0
0.0   1.0 5.0
0.0   2.0 5.0
0.0   3.0 5.0
0.0   4.0 5.0
0.0   5.0 5.0
0.0   6.0 5.0
0.0   7.0 5.0
0.0   8.0 5.0
0.0   9.0 5.0
0.0  10.0 5.0
plume line0
Figure 3. An extended plume (in Blender) from a straight line of emitters with equal probabilities.

The example file below shows the same set of 21 point emitters as the previous example, but this example has a probability associated with each emitter. If a 4th token is not present in a given emitter description then the probability is assumed to be 1. This example has the probabilities lower at the ends of the extended source with a lower probability near the center as well.

An example extended source emitter file containing 21 locations along a straight line with varying probabilities.
0.0 -10.0 5.0 0.2
0.0  -9.0 5.0 0.3
0.0  -8.0 5.0 0.4
0.0  -7.0 5.0 0.5
0.0  -6.0 5.0 0.8
0.0  -5.0 5.0 1.0
0.0  -4.0 5.0 1.0
0.0  -3.0 5.0 1.0
0.0  -2.0 5.0 1.0
0.0  -1.0 5.0 0.4
0.0   0.0 5.0 0.1
0.0   1.0 5.0 0.4
0.0   2.0 5.0 1.0
0.0   3.0 5.0 1.0
0.0   4.0 5.0 1.0
0.0   5.0 5.0 1.0
0.0   6.0 5.0 0.8
0.0   7.0 5.0 0.5
0.0   8.0 5.0 1.4
0.0   9.0 5.0 0.3
0.0  10.0 5.0 0.2
plume line1
Figure 4. An extended plume (in Blender) from a straight line of emitters with varying probabilities.

The example file below shows a set of 21 point emitters that are primarily distributed along the Y-axis, with some variation in the X-axis position and a constant Z-axis location.

An example extended source emitter file containing 21 locations along an irregular line with probabilities.
 0.0 -10.0 5.0 0.2
-1.0  -9.0 5.0 0.3
-1.0  -8.0 5.0 0.4
-2.0  -7.0 5.0 0.5
-4.0  -6.0 5.0 0.8
-3.0  -5.0 5.0 1.0
-7.0  -4.0 5.0 1.0
-2.0  -3.0 5.0 1.0
-1.0  -2.0 5.0 1.0
 0.0  -1.0 5.0 0.4
 0.0   0.0 5.0 0.1
 0.0   1.0 5.0 0.4
-1.0   2.0 5.0 1.0
-2.0   3.0 5.0 1.0
-3.0   4.0 5.0 1.0
-5.0   5.0 5.0 1.0
-5.0   6.0 5.0 0.8
-3.0   7.0 5.0 0.5
-1.0   8.0 5.0 1.4
-1.0   9.0 5.0 0.3
 0.0  10.0 5.0 0.2
plume line2
Figure 5. An extended plume (in Blender) from an irregular line of emitters with equal probabilities.

Notes

  • With DIRSIG5, the plugin instantiates the plume into the scene so it is recommended 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 0.1 second intervals. the stepcount can be used to compute 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 = 0.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 artificial 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 tool in the DIRSIG vdb_tool utility program is run by supplying the name of the input XML file and the name of the output VDB file:

Example usage of the blackadar tool.
$ vdb_tool blackadar -o plume.vdb plume.xml

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 tool with a random seed.
$ vdb_tool blackadar --seed=1234 --output=plume.vdb plume.xml

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.