Lab 3-3: Stability classes#
# to try this out, let's return to our Gaussian plume model but with different stabilities
# first, load our packages and read in our data.
import numpy as np
import scipy.stats as st
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
# In this lab, we will be using if and elif statements to assign variables.
test = 'AA'
if test == 'AA':
print('This worked')
This worked
# Let's start with the functions and parameters from the Martin 1976 paper
# here, we focus on distances more than 1 km (1000 m) from the source
# for shorter distances, there are different functional forms.
# Note that all of these are for vertical spreading sigma_z -- we will focus
# our analysis along the centerline y (although the paper has values for sigma_y as well)
# specify the range of interest over x (>1000 m only)
x = np.arange(1000,100000,1)
# b is a constant value across stability classes
b = 0.894
# specify the atmospheric class (based on lecture notes and Martin 1976 paper)
airclass = 'A'
if airclass == 'A':
c = 459.7
d = 2.094
f = -9.6
sig_z = c*np.power(x/1000,d) + f
a = 213
sig_y = a*np.power(x/1000,b)
elif airclass == 'B':
c = 108
d = 1.1
f = 2.0
sig_z = c*np.power(x/1000,d) + f
a = 156
sig_y = a*np.power(x/1000,b)
elif airclass == 'C':
c = 61
d = 0.911
f = 0
sig_z = c*np.power(x/1000,d) + f
a = 104
sig_y = a*np.power(x/1000,b)
elif airclass == 'D':
c = 44.5
d = 0.516
f = -13
sig_z = c*np.power(x/1000,d) + f
a = 68
sig_y = a*np.power(x/1000,b)
elif airclass == 'E':
c = 55.4
d = 0.305
f = -34
sig_z = c*np.power(x/1000,d) + f
a = 50.5
sig_y = a*np.power(x/1000,b)
elif airclass == 'F':
c = 62.6
d = 0.180
f = -48.6
sig_z = c*np.power(x/1000,d) + f
a = 34
sig_y = a*np.power(x/1000,b)
# print out the value of c to make sure that the program assigned the right one
print(c)
459.7
#Let's try plotting how sigma changes with distance downwind, x
plt.figure()
plt.plot(x, sig_z, label=airclass)
plt.ylabel('sigma')
plt.xlabel('distance downwind (m)')
plt.title('Smokestack spreading function')
plt.legend(loc="best")
<matplotlib.legend.Legend at 0x7f33751a7390>
# Following the lecture notes, we can then plot the ground concentration along the centerline
# where y = 0 and z = 0.
# if we assume the same plume discharge, height and wind as in lab 3-1, we get:
Q = 8500
# Q is a constant mass flux from the smokestack in kg/s
u = 3.5
# u is a constant wind speed in m/s at the effective stack height h
h = 200
# h is the effective height of the smokestack (actually stack height plus additional
# height the plume rises, in m
C=Q/(np.pi*u*sig_z*sig_y)*np.exp(-(np.power(h,2))/(2*np.power(sig_z,2)))
plt.figure()
plt.plot(x, C, label=airclass)
plt.ylabel('kg/m^3')
plt.xlabel('distance downwind (m)')
plt.title('Plume concentration at centerline along the ground')
plt.legend(loc="best")
<matplotlib.legend.Legend at 0x7f3374f5d190>
# We can also make a plot where we zoom in closer (as concentration is very negligible far downwind
plt.figure()
plt.plot(x, C, label=airclass)
plt.ylabel('kg/m^3')
plt.xlabel('distance downwind (m)')
plt.xlim(0,10000)
plt.title('Plume concentration at centerline along the ground')
plt.legend(loc="best")
<matplotlib.legend.Legend at 0x7f3375034d10>
Work on your own#
Now, try changing the stability class and comparing the plots. You may want to assign different variable names to the different spreading functions and concentrations so that you can save them all and then plot them all on one graph for comparison.
# Let's plot the plume concentration for three different stability classes to compare.
# start with presuming the same wind speed in all cases for comparison
u100_1= 5 # Case 1, F, very stable
u100_2= 5 # Case 2, D, neutral
u100_3= 5 # Case 3, A, unstable
# if we assume the same plume discharge, height and wind as in lab 3-1, we get:
Q = 8500
# Q is a constant mass flux from the smokestack in kg/s
h = 100
# h is the effective height of the smokestack (actually stack height plus additional
# height the plume rises, in m
b = 0.894 # Martin 1976 defines this as the same across stability classes
#airclass == 'F':
c = 62.6
d = 0.180
f = -48.6
sig_z = c*np.power(x/1000,d) + f
a = 34
sig_y = a*np.power(x/1000,b)
u = u100_1
# u is a constant wind speed in m/s at the effective stack height h, and will vary with stability class
C1=Q/(np.pi*u*sig_z*sig_y)*np.exp(-(np.power(h,2))/(2*np.power(sig_z,2)))
#airclass == 'D':
c = 44.5
d = 0.516
f = -13
sig_z = c*np.power(x/1000,d) + f
a = 68
sig_y = a*np.power(x/1000,b)
u = u100_2
# u is a constant wind speed in m/s at the effective stack height h, and will vary with stability class
C2=Q/(np.pi*u*sig_z*sig_y)*np.exp(-(np.power(h,2))/(2*np.power(sig_z,2)))
#airclass == 'A':
c = 459.7
d = 2.094
f = -9.6
sig_z = c*np.power(x/1000,d) + f
a = 213
sig_y = a*np.power(x/1000,b)
u = u100_3
# u is a constant wind speed in m/s at the effective stack height h, and will vary with stability class
C3=Q/(np.pi*u*sig_z*sig_y)*np.exp(-(np.power(h,2))/(2*np.power(sig_z,2)))
plt.figure()
plt.plot(x, C1, label='Case 1')
plt.plot(x, C2, label='Case 2')
plt.plot(x, C3, label='Case 3')
plt.ylabel('kg/m^3')
plt.xlabel('distance downwind (m)')
plt.xlim(0,10000)
plt.title('Plume concentration at centerline along the ground')
plt.legend(loc="best")
<matplotlib.legend.Legend at 0x7f3374cb7250>
# and now I want to add the vertical component, including the reflection
# consider the centerline
y=0
u=5
# and compare plots for 1000 m for the classes, and then, a different plot for 2000 m downwind
z = np.arange(0,1000,1)
dist_km=1 #change to a different distance downwind to plot
#airclass == 'F':
c = 62.6
d = 0.180
f = -48.6
sig_z = c*np.power(dist_km,d) + f
a = 34
sig_y = a*np.power(dist_km,b)
Cwreflect1=Q/(2*np.pi*u*sig_z*sig_y)*(np.exp(-(np.power(z-h,2))/(2*np.power(sig_z,2)))+np.exp(-(np.power(z+h,2))/(2*np.power(sig_z,2))))*np.exp(-(np.power(y,2))/(2*np.power(sig_y,2)))
#airclass == 'D':
c = 44.5
d = 0.516
f = -13
sig_z = c*np.power(dist_km,d) + f
a = 68
sig_y = a*np.power(dist_km,b)
Cwreflect2=Q/(2*np.pi*u*sig_z*sig_y)*(np.exp(-(np.power(z-h,2))/(2*np.power(sig_z,2)))+np.exp(-(np.power(z+h,2))/(2*np.power(sig_z,2))))*np.exp(-(np.power(y,2))/(2*np.power(sig_y,2)))
#airclass == 'A':
c = 459.7
d = 2.094
f = -9.6
sig_z = c*np.power(dist_km,d) + f
a = 213
sig_y = a*np.power(dist_km,b)
Cwreflect3=Q/(2*np.pi*u*sig_z*sig_y)*(np.exp(-(np.power(z-h,2))/(2*np.power(sig_z,2)))+np.exp(-(np.power(z+h,2))/(2*np.power(sig_z,2))))*np.exp(-(np.power(y,2))/(2*np.power(sig_y,2)))
plt.figure()
plt.plot(Cwreflect1, z, label='Case 1')
plt.plot(Cwreflect2, z, label='Case 2')
plt.plot(Cwreflect3, z, label='Case 3')
plt.xlabel('kg/m^3')
plt.ylabel('elevation (m)')
plt.ylim(0,300)
plt.title('Plume concentration with height at 1 km downwind')
plt.legend(loc="best")
<matplotlib.legend.Legend at 0x7f3374d1a110>
# and now I want to add the vertical component, including the reflection
# consider the centerline
y=0
u=5
# and compare plots for 1000 m for the classes, and then, a different plot for 2000 m downwind
z = np.arange(0,1000,1)
dist_km=2 #change to a different distance downwind to plot
#airclass == 'F':
c = 62.6
d = 0.180
f = -48.6
sig_z = c*np.power(dist_km,d) + f
a = 34
sig_y = a*np.power(dist_km,b)
Cwreflect1=Q/(2*np.pi*u*sig_z*sig_y)*(np.exp(-(np.power(z-h,2))/(2*np.power(sig_z,2)))+np.exp(-(np.power(z+h,2))/(2*np.power(sig_z,2))))*np.exp(-(np.power(y,2))/(2*np.power(sig_y,2)))
#airclass == 'D':
c = 44.5
d = 0.516
f = -13
sig_z = c*np.power(dist_km,d) + f
a = 68
sig_y = a*np.power(dist_km,b)
Cwreflect2=Q/(2*np.pi*u*sig_z*sig_y)*(np.exp(-(np.power(z-h,2))/(2*np.power(sig_z,2)))+np.exp(-(np.power(z+h,2))/(2*np.power(sig_z,2))))*np.exp(-(np.power(y,2))/(2*np.power(sig_y,2)))
#airclass == 'A':
c = 459.7
d = 2.094
f = -9.6
sig_z = c*np.power(dist_km,d) + f
a = 213
sig_y = a*np.power(dist_km,b)
Cwreflect3=Q/(2*np.pi*u*sig_z*sig_y)*(np.exp(-(np.power(z-h,2))/(2*np.power(sig_z,2)))+np.exp(-(np.power(z+h,2))/(2*np.power(sig_z,2))))*np.exp(-(np.power(y,2))/(2*np.power(sig_y,2)))
plt.figure()
plt.plot(Cwreflect1, z, label='Case 1')
plt.plot(Cwreflect2, z, label='Case 2')
plt.plot(Cwreflect3, z, label='Case 3')
plt.xlabel('kg/m^3')
plt.ylabel('elevation (m)')
plt.ylim(0,300)
plt.title('Plume concentration with height at 2 km downwind')
plt.legend(loc="best")
<matplotlib.legend.Legend at 0x7f3374bacd10>
