In [1]:
from IPython.display import Image

Weighting functions in the $CO_2$ 15 $\mu m$ absorption band

Below is a plot of radiance (or intensity) (left axis) and brightness temperature (right axis) vs. wavenumber near the main $CO_2$ absorption band. Wavenumber is defined as $1/\lambda$; the center of the absorption band is at $\lambda = 15\ \mu m$ which is a wavenumber of 1/0.0015 = 666 $cm^{-1}$. The VTPR (vertical temperature profiling radiometer) has six channels: Channel 1 is at the center of the band -- it has the lowest transmissivity and is measuring photons coming from 45 km at the top of the stratosphere (See Stull Chapter 1, Figure 1.10). As the channel number increases from 2-6 the transmissivity also increases, and the photons originate from increasing lower levels of the atmosphere with increasing kinetic temperatures. Note that the different heights for the peaks in each of the weighting functions.


In [2]:
Image('figures/wallace4_33.png',width=500)


Out[2]:

Assignment 9

In this assignmet we'll work with the five standard atmospheres from hydrostatic.ipynb to show how to use Stull eq. 8.4 to calculate radiance at the top of the atmosphere for wavelengths similar to the 6 sounder channels shown in the above figure. I define a new function find_tau below to calculate the optical thickness for a $CO_2$-like absorbing gas. I ask your to find the transmissivities $t=\exp(-\tau_{tot} - \tau)$ and the weighting functions $\Delta t$, at 7 wavelengths, and use those plus $B_\lambda$ to calculate the radiance for the 7 channels for each of the 5 atmospheres.


In [3]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticks
from a301utils.a301_readfile import download
import h5py
from a301lib.radiation import Blambda,planckInvert
#
# from the hydrostatic noteboo
#
from a301lib.thermo import calcDensHeight

In [4]:
filename='std_soundings.h5'
download(filename)


std_soundings.h5 already exists, need to either delete or rename for new download

In [5]:
from pandas import DataFrame
with h5py.File(filename) as infile:
    sound_dict={}
    print('soundings: ',list(infile.keys()))
    #
    # names are separated by commas, so split them up
    # and strip leading blanks
    #
    column_names=infile.attrs['variable_names'].split(',')
    column_names = [item.strip() for item in column_names]
    column_units = infile.attrs['units'].split(',')
    column_units = [item.strip() for item in column_units]
    for name in infile.keys():
        data = infile[name][...]
        sound_dict[name]=DataFrame(data,columns=column_names)


soundings:  ['midsummer', 'midwinter', 'subsummer', 'subwinter', 'tropics']

mass absorption coefficient for fake $CO_2$

To keep things simple I'm going to make up a set of 7 absorption coefficients that will give weighting functions that look something like the VPTR.

We have been working with the volume absorption coefficient:

$$\beta_\lambda = n b\ (m^{-1})$$

where $n$ is the absorber number concentration in $\#\,m^{-3}$ and $b$ is the absorption cross section of a molecule in $m^2$. For absorbing gasses like $CO_2$ that obey the ideal gas law $n$ depends inversely on temperature -- which changes rapidly with height.

A better route is to use the mass absorption coefficient $k_\lambda$:

$$k_\lambda = \frac{n b}{\rho_{air}}$$

units: $m^2/kg$. With this definition the optical depth is:

$$\tau_\lambda = \int_0^z \rho_{air} k_\lambda dz$$

Why is this an improvement? We now have an absorption coefficient that can is roughly constant with height and can be taken out of the integral.


In [6]:
Rd = 287.
def find_tau(r_gas,k_lambda,df):
    """
    given a data frame df with a standard sounding, return the optical depth assuming an
    absorbing gas with a constant mixing ration r_gas and absorption coefficient k_lambda
    
    Parameters
    ----------
    
    r_gas: float
       absorber mixing ratio, kg/kg
       
    k_lambda: float
       mass absorption coefficient, m^2/kg
       
    df:  dataframe
       sounding with height levels as rows, columns 'temp': temperature in K
       'press': pressure in Pa
       
    Returns
    -------
    tau: vector (float)
      optical depth measured from the surface at each height in df
    """
    #
    # density scale height
    #
    Hdens = calcDensHeight(df)
    #
    # surface density
    #
    rho0 = df['press']/(Rd*df['temp'])
    rho = rho0*np.exp(-df['z']/Hdens)
    height = df['z'].values
    tau=np.empty_like(rho)
    #
    # start from the surface
    #
    tau[0]=0
    num_levels=len(rho)
    num_layers=num_levels-1
    for index in range(num_layers):
        delta_z=height[index+1] - height[index]
        delta_tau=r_gas*rho[index]*k_lambda*delta_z
        tau[index+1]=tau[index] + delta_tau
    return tau

Example $\tau$ calculation for r_gas=0.01 $kg/kg$ and k_lambda = 0.01 $m^2/kg$


In [7]:
%matplotlib inline
r_gas=0.01  #kg/kg
k_lambda=0.01  #m^2/kg
df=sound_dict['tropics']
#
# set the top of the atmosphere at 20 km
#
top = 20.e3
df = df.loc[df['z']<top]
height = df['z']
press = df['press']
tau=find_tau(r_gas,k_lambda,df)
fig1,axis1=plt.subplots(1,1)
axis1.plot(tau,height*1.e-3)
axis1.set_title('vertical optical depth vs. height')
axis1.set_ylabel('height (km)')
axis1.set_xlabel('optical depth (no units)')
fig2,axis2=plt.subplots(1,1)
axis2.plot(tau,press*1.e-3)
axis2.invert_yaxis()
axis2.set_title('vertical optical depth vs. pressure')
axis2.set_ylabel('pressure (kPa)')
axis2.set_xlabel('optical depth (no units)')


Out[7]:
<matplotlib.text.Text at 0x115bb9ac8>

Extending this to 7 wavelengths

In the next cell I make up 7 k_lambda values to go with 7 wavelengths from 13 to 15 microns (766 to 666 $cm^{-1}$)


In [8]:
#
# assign the 7 k_lambdas to 7 CO2 absorption band wavelengths
# (see Wallace and Hobbs figure 4.33)
#
wavenums=np.linspace(666,766,7)  #wavenumbers in cm^{-1}
wavelengths=1/wavenums*1.e-2  #wavelength in m
wavelengths_um = wavelengths*1.e6  # in microns
print('channel wavelengths (microns) ',wavelengths_um)  #microns
df=sound_dict['tropics']
top = 20.e3  #stop at 20 km
df = df.loc[df['z']< top]
height = df['z'].values
mid_height = (height[1:] - height[:-1])/2.
#
# here are the mass absorption coefficients for each of the 7 wavelengths
# in m^2/kg
#
k_lambda_list=np.array([ 0.175 , 0.15 ,  0.125 , 0.1  ,  0.075,  0.05 ,  0.025])
legend_string=["{:5.3f}".format(item) for item in k_lambda_list]
#
#  find the height at mid-layer
#
mid_height=(height[1:] + height[:-1])/2.
#
# make a list of tuples of k_lambda and its label
# using zip
#
k_vals=zip(k_lambda_list,legend_string)
fig1,ax=plt.subplots(1,1,figsize=(10,10))
heightkm=height*1.e-3
mid_heightkm=mid_height*1.e-3
for k_lambda,k_label in k_vals:
    tau=find_tau(r_gas,k_lambda,df)
    ax.plot(tau,heightkm,label=k_label)
    ax.legend()


channel wavelengths (microns)  [ 15.01501502  14.6484375   14.2993327   13.96648045  13.64877161
  13.34519573  13.05483029]

To hand in

1) make two more plots for the tropical atmosphere:

a) a plot of the transmissivity $t_\lambda = exp(-(\tau_{tot} - \tau)$ for each wavelength vs. height in km, with a legend, axes labels and a title

b) a plot of the weighting functions $\Delta t_\lambda$ vs. height for each wavelength, with a legend, axes labels and a title

2) Use Stull 8.4 with Blambda and planckInvert imported above to calculate the brighntess temperature for each of the 5 atmospheres at each of the 7 wavelengths. Make a plot of wavelength vs. brightness temperature with 5 lines on it, one for each atmosphere. The plot should have a title, axes labels and a legend for the 5 lines

Hint: I wrote two nested loops, looping through the dataframes in sound_dict and then looping over the 7 wavelengths with 7 k_lambdas giving me 7 taus, transmissivities and weighting functions. I initialized a dictionary with each sounding name as a key containing empty lists to hold the 7 brightness temperatures. Then in my wavelength loop I used Stull 8.4 to calculate the radiance at 20 km (our top of atmosphere), inverted to get brightness temperature and appended that to the list for that sounding. The brightness temperatures should range between about 220 K and 290 K


In [ ]: