This summer has been a bad one for wildfires and the resultant smoke filling the city, today is no exception and the city is blanketed in a thick cloud of smoke. Thick enough that, standing on my balcony, I can barely see buildings in my neighbourhood through the gloom. I've taken to habitually checking the smoke forecast before going outside and doing anything. In the process it got me thinking about dispersion models and how they work. Beyond forecasting the forest fire smoke, they are used throughout industry to predict where emissions from vents and stacks will go, and also in modeling emergency scenarios (if there is an uncontrolled release of some kind, where will it go?).

So let's imagine that I'm standing next to a burning tree with an equally burning desire to know where the smoke will go. There are many complex models that take into account the surface features of the earth and the dynamics of the atmosphere in calculating where that smoke will go, but I'm just sitting in the woods with a notebook and a calculator, so where to start?

Suppose we start with some differential volume element in the air, we can do a mass balance over the element fairly easily:

$$ Accumulation = In - Out $$ $$ {\partial C \over \partial t} V = \left( N_{x} A_{x} - N_{x+dx} A_{x+dx} \right) + \left( N_{y} A_{y} - N_{y+dy} A_{y+dy} \right) + \left( N_{z} A_{z} - N_{z+dz} A_{z+dz} \right)$$

Where *C* is the mass concentration of smoke particles, $N_{x}$ is the mass flux through the surface at *x*, and $N_{x+dx}$ is the mass flux through the surface at *x + dx*.

The total mass in at *x* is a sum of diffusion (Fick's Law) and bulk flow contributions.

$$ N_{x} \> dy dz = - {\partial \over \partial x} D_{x} C \> dy dz + u_{x} C \> dy dz$$

Where *D* the diffusivity and *u* the fluid velocity.

Since this is a small (differential) element we can linearly approximate the total mass flow at *x + dx*:

$$ N_{x + dx} \> dy dz = N_{x} \> dy dz + {\partial \over \partial x} \left ( N_{x} \> dy dz \right ) \> dx $$ $$ N_{x + dx} \> dy dz = - {\partial \over \partial x} D_{x} C \> dy dz + u_{x} C \> dy dz + {\partial \over \partial x} \left ( - {\partial \over \partial x} D_{x} C \> dy dz + u_{x} C \> dy dz \right ) \> dx $$

Similarly for y and z.

Let's assume we've oriented our coordinate system such that the air is moving entirely in the x direction, i.e. $u_{y} = u_{z} = 0$ we can put everything together, do some simplifying and we get:

$$ {\partial C \over \partial t} \> dx dy dz = \left( - {\partial \over \partial x} D_{x} C \> dy dz + u_{x} C \> dy dz + {\partial \over \partial x} D_{x} C \> dy dz - u_{x} C \> dy dz - {\partial \over \partial x} \left ( - {\partial \over \partial x} D_{x} C \> dy dz + u_{x} C \> dy dz \right ) \> dx \right) + \ldots $$

$$ {\partial C \over \partial t} \> dx dy dz = \left( {\partial^{2} \over \partial x^{2}} D_{x} C - {\partial \over \partial x} u_{x} C \right ) \> dx dy dz + \left( {\partial^{2} \over \partial y^{2}} D_{y} C \right ) \> dx dy dz + \left( {\partial^{2} \over \partial z^{2}} D_{z} C \right ) \> dx dy dz$$

$$ {\partial C \over \partial t} = - {\partial \over \partial x} u_{x} C + \left( {\partial^{2} \over \partial x^{2}} D_{x} C + {\partial^{2} \over \partial y^{2}} D_{y} C + {\partial^{2} \over \partial z^{2}} D_{z} C \right )$$

Now for some more simplifying assumptions:

- Assume (with some justification) that $ \left | {\partial^{2} \over \partial x^{2}} D_{x} C \right | \ll \left | {\partial \over \partial x} C u \right | $ basically that the wind is much more significant than diffusion in the x direction.
- The system is at steady state, so ${\partial C \over \partial t} = 0$
- The wind speed
*u*is spatially constant (this is very untrue) - The
*D*'s are spatially constant

Throwing this all into a blender we end up with:

$$ u {\partial C \over \partial x} = D_{y} {\partial^{2} C \over \partial y^{2}} + D_{z} {\partial^{2} C \over \partial z^{2}} $$

Which is decidedly more solvable, in fact we can just look up the answer (after having made some reasonable boundary condition assumptions):

$$ C = {k \over x} \exp \left[ - \left( {y^{2} \over D_{y}} + {z^{2} \over D_{z}} \right) { u \over 4x } \right] $$

Where *k* is a constant set by the boundary conditions.

At this point we take a deep breath and look up from our notebook, is the tree still burning? Is the whole forest burning? Is this mathematical interlude madness and folly? What have we ended up with? We have a model of a steady state plume of smoke, that starts at a source and is swept downwind, dispersing outward as we go. We haven't accounted for any flow due to buoyancy or anything else. It's very simple, but hopefully representative.

Don't get discouraged! We only need to solve for *k*!

Ok so what do we know? We've assumed that the system is at a steady state and that there is no diffusion in the x direction (merely bulk flow). If we suppose a constant amount of smoke is produced, *Q* kg/s of "smoke", then that mass has to be conserved as it is carried downwind. If we imagine slicing our smoke plume into slices along the yz plane then the mass that passes through each slice per second is *Q* (and this is the same for any slice).

$$ Q = \int \int C u \> dy dz $$ $$ Q = \int_{0}^{\infty} \int_{-\infty}^{\infty} {k u \over x} \exp \left[ - \left( {y^{2} \over D_{y}} + {z^{2} \over D_{z}} \right) { u \over 4x } \right] \> dy dz $$

Where we imagine (for simplicity) that the fire is at ground level (z = 0).

Let $y' = {y \over \sqrt{D_{y}}}$ and $z' = {z \over \sqrt{D_{z}}}$ then:

$$ Q = {k u \over x} \sqrt{D_{y} D_{z}} \> \int_{-\infty}^{\infty} \exp \left[ - {u \over 4 x} y'^{2} \right] \> dy' \> \int_{0}^{\infty} \exp \left[ - {u \over 4 x} z'^{2} \right] \> dz' $$

These are now two gaussian integrals, which we can look up, to get:

$$ Q = {k u \over x} \sqrt{D_{y} D_{z}} \> \left( \sqrt{\pi} \over \sqrt{u \over 4x} \right) \> \left( \sqrt{\pi} \over 2 \sqrt{u \over 4x} \right) $$ $$ Q = 2 \pi k \sqrt{D_{y} D_{z}} $$

$$k = {Q \over 2 \pi \sqrt{D_{y} D_{z}} }$$

Deep breath! Finally we have:

$$ C = {Q \over 2 \pi x \sqrt{D_{y} D_{z}} } \exp \left[ - \left( {y^{2} \over D_{y}} + {z^{2} \over D_{z}} \right) { u \over 4x } \right] $$

What we just derived is a *Gaussian Dispersion Model*. It is the first model you try when modeling releases in the chemical industry, from smoke stacks or other point sources. In this case I assumed the point source was on the ground for simplicity, in more developed models the point source is at some height *H* and an effect called "ground reflection" needs to be taken into account.

An alternative way of writing the model is to let $ \sigma_{y}^{2} = {2 D_{y} x \over u}$ and $ \sigma_{z}^{2} = {2 D_{z} x \over u}$ (making it more explicitly *gaussian* in form) and write the equation:

$$ C = {Q \over \pi u \sigma_{y} \sigma_{z} } \exp \left[ -\frac{1}{2} \left( {y^{2} \over \sigma_{y}^{2}} + {z^{2} \over \sigma_{z}^{2}} \right) \right] $$

Where $\sigma_{y}$ and $\sigma_{z}$ are either tabulated, as a functon of atmospheric stability, or provided through some empirical correlation. My textbook gives a functional form:

$$ \sigma_{y} = a x^{b} $$ $$ \sigma_{z} = c x^{d} + f $$

With constants tabulated based on the Pasquill stability class criteria, and x is the distance downwind of the point source.

One big limitation of this model is that it assumes the wind is contant, over time and over space, which not true in general. So this gives a nice *impression* of what the concentration will be downwind of the point source, but in general it will only give an impression. To go further one must account for the dynamic nature of the environment (i.e. the weather) which multiplies the complexity a lot!