In the last article we looked at some results from Grant Petty.
In essence they demonstrate the huge variation in the transmittance of CO2 through the atmosphere with the wavelength of radiation.
I decided it might be interesting to try and reproduce these results using the HITRAN database. This allows a closer examination of which mechanisms cause which results.
Line Shapes and Pressure Broadening
The energy in a photon is given by a very simple equation:
E = hν
where h = Planck’s constant = 6.63 x 10-34 J.s, and ν = frequency (ν=c/λ where c=speed of light and λ = wavelength)
So, for example, a 15μm photon has an energy of 1.3 x 10-20 J.
A photon can be absorbed by a molecule if the energy of the photon matches the exact energy required to change the state of that molecule.
Now, if a photon was only absorbed at one exact wavelength then atmospheric radiation would be irrelevant.
An easy concept for people who have done a lot of maths or physics, but not so obvious to many others.
Let’s take a different example. If you start your car up and accelerate steadily up to 60 miles per hour over 1 minute, how long do you spend at the speed of 34.5698895549034592345123 miles per hour? Not much. And the more precise I make this speed, the less time you will have spent at it.
If we consider how much energy is transferred in terrestrial radiation from the surface at 14.995698895549034592345123 μm, the answer is “not much”. The same goes for any other vanishingly small interval of wavelengths.
So if each absorption line was exactly one wavelength then the amount of energy absorbed by a finite number of lines would be zero.
However, each absorption line has a finite width. Part of this is due to the uncertainty principle described by quantum mechanics – the small time a higher energy state is occupied produces an uncertainty in the energy of that state.
This is called natural broadening and is a very small effect.
The easiest broadening mechanism to understand is Doppler broadening. Molecules in the atmosphere are zipping around in all directions at speeds of around 500 m/s. Compared with the speed of light this is small, but the difference in the relative speed between a photon and a molecule is slightly changed, leading to a change in the absorption frequency (see note 1 for more details). This is “relatively easy” to understand from first principles and the mathematical treatment is straightforward.
The dominant line broadening mechanism in the lower atmosphere is pressure, or collisional, broadening. This is more challenging from a theoretical point of view but is certainly a measurable value.
The main effect of pressure broadening is to “spread out” the absorption line. Here is a typical example, where 1000mb is at the earth’s surface, and 200mb is at the “tropopause”, or top of the lower atmosphere:
So the absorption at the center of the line is reduced, while more absorption takes place out to the sides.
The parameter which describes how much broadening has taken place is called the half-width – and is the width of the line when the value has dropped to half the peak value.
Typical values of pressure broadening half width are 0.01 – 0.1 cm-1 at 296 K and pressure of 1000 mb (surface atmospheric pressure).
The function which approximates this effect is the Lorentz line shape:
f(ν-ν0) = αL/π((ν-ν0)²+αL²) 
where ν0 is the frequency of the absorption, ν = frequency and αL = half-width at half of the maximum
And the half-width at any given temperature and pressure is given by this formula:
αL = α0(p/p0).(T0/T)γ 
where α0 = lab measured half-width at pressure p0 and temperature T0, and γ is also a lab measured parameter, typically around 0.5 – 0.7.
That’s a lot to take in if you haven’t seen it before. Take a look back at Figure 1 – this is what these formulas describe.
The atmospheric pressure where pressure broadening is comparable to Doppler broadening is about 10 mbar. This is around 30km above the surface, so in the troposphere Doppler broadening is unimportant.
In a later article we might explore the experimental results vs theoretical results of pressure broadening in more detail. Or at least, how much inaccuracies might affect radiative transfer calculations.
You can read more about the HITRAN database here. The most recent update of the database is described in The HITRAN 2008 molecular spectroscopic database, by L.S. Rothman et al, and the earlier 2004 update: The HITRAN 2004 molecular spectroscopic database, by L.S. Rothman et al.
This database is the result of a huge amount of work by thousands of researchers over a few decades. If you look at the references in the 2008 paper you find almost 400 papers.
In the previous articles in this series I created two fictitious molecules, pCO2 and pH20, and solved the radiative transfer equations for a variety of conditions for these two molecules through the atmosphere.
The aim of this simplification was illumination.
Given that the HITRAN database contains over 300,000 absorption lines for CO2 (and 2.7M lines in total), I thought that replicating some standard results (e.g. see Part Eight) might be a little too much of a challenge. However, the hardest part was actually reading the database, which is in text form – and the problems were just due to my novice status with MATLAB.
DeWitt Payne kindly offered me some results he obtained from SpectralCalc (via subscription). These were from 660-672 cm-1 through a 1 meter path at standard temperature and pressure, and CO2 concentrations of 380 ppm.
Here is the comparison with my MATLAB program:
And here is the difference between the two results (the SpectralCalc results minus the MATLAB results):
The MATLAB program, for the results above, only considered the main isotopologue of CO2, which accounts for over 98% of the concentration of CO2 in the atmosphere. However, including the absorption lines for all isotopologues had almost no effect and the small differences remain.
As the differences are small, it seems worth showing some initial results from the program.
First, let’s see the effect of “pressure broadening”.
Here is the result through 1m of atmosphere at the earth’s surface, around the peak absorption wavenumber of CO2:
And the result for the same 1m of atmosphere at the tropopause (top of the troposphere) at 200mbar (typically around 12km):
Now at 200 mbar there is less atmosphere so the comparison isn’t truly a fair one. Increasing the path length to get the same number of CO2 molecules we get this result:
Both pressure and temperature have changed so it’s worth seeing the result if we apply the surface temperature to the Figure 5 results:
We can see from the result (and it’s also clear from the equations earlier) that the parameter having the most effect is the pressure – which changes by a factor of 5. The temperature reduction from 296K to 216K has a much smaller impact.
However, the difference between Figure 4 and Figure 6 is very significant – individual lines are “smeared out” at the earth’s surface, but distinct lines higher up in the atmosphere.
Let’s take a look at a wider range of wavenumbers / wavelengths.
From 600 – 750 cm-1 (13.3 – 16.7 μm) at the surface through 1 meter of atmosphere:
And the same amount of CO2 at the tropopause:
Let’s look at a longer path of 100m through the atmosphere and “zero in” on one part of the bandwidth, 640-650 cm-1:
And let’s compare it with the tropopause, again through a longer path, to keep the CO2 quantity the same:
Just for interest I compared the optical thickness of the surface and tropopause for these wavenumbers on the same graph. Optical thickness increases as transmittance decreases. Check out the heading Optical Thickness & Transmittance in Part Six if this isn’t clear.
We can see that the peaks of absorption are indeed lower but the widths of the absorption lines are wider (for the surface results).
As I commented in Part Eight, many people have heard about the high absorption of CO2 at 15 μm and have not understood the huge variation in absorption across the wide bandwidth of CO2. (Also there is the most important subject of re-emission).
Here is the result of the simple “slab” model, which calculates the transmittance through 1km of atmosphere (a slab model means that the pressure and temperature are constant through the “slab” – a very simple model):
As you can see, around 570 – 600 cm-1 (16.7 – 17.5 μm) and 730 – 770 cm-1 (13.0 – 13.7 μm) the transmittance through the atmosphere is nowhere near “saturated”.
If 1 km of atmosphere has a transmittance of 0.7 (i.e., 70% of radiation is transmitted) at some wavelengths then clearly more CO2 will reduce this transmittance.
The Matlab Code
The results were all created using my code: HITRAN_0_2. However, each run has slightly different parameters, which are adjusted by changing the code rather than anything more elegant. See Note 2 for the code.
The code at this stage has provision for different layers of different temperatures, pressures, and densities – but at this stage the code doesn’t use it – just one “slab”.
The core part of the code is very simple:
- Take each absorption line in turn and use its measured line strength, half width, and other parameters
- For the pressure and temperature, calculate the half width for those conditions (equation 2)
- Now for each line, take each wavenumber in turn (defined via a start, stop and interval) and apply the formula in equation 1 to calculate the optical thickness
- Add the optical thickness at this wavenumber to the optical thickness already calculated for this wavenumber
The code could be more “speed efficient” – by only calculating equation 1 close to the absorption line. I started out simpler to see what happened and it went so fast that I didn’t improve it.
Multiple Atmospheric Layers
As the pressure and temperature change with height through the atmosphere the single slab model is limited in value.
So I extended the model to multiple layers – v0.3 (see Note 2).
Here is the result of transmittance up to 200mb (20,000 Pa) from CO2 at 360 ppm with a 15 layer model:
And, of course, the question many people have – what is the difference between the atmosphere at 280 ppm and 560 ppm – or a doubling of CO2 from pre-industrial levels:
And the difference between the two, the “delta”:
Remember – or note, if this is new – the atmosphere absorbs and also re-emits. The question of radiative forcing cannot be answered by only considering transmittance.
However, these graphs of the change should at least indicate to people that the issue of the very high optical thickness = very low transmittance at 667 cm-1 = 15 μm is not the complete story, even only considering transmittance.
The model has the limitation at the moment that it doesn’t take into account Doppler broadening. If we extend the model up into the stratosphere Doppler broadening starts to become important. I might get around to including this by introducing the Voigt function which combines pressure broadening and Doppler broadening.
There is also the possibility that my model has a mistake. However, I have compared it to the results from SpectralCalc and it is very close, and to results published in Grant Petty’s book, which look the same.
Part One – a bit of a re-introduction to the subject.
Part Two – introducing a simple model, with molecules pH2O and pCO2 to demonstrate some basic effects in the atmosphere. This part – absorption only.
Part Three – the simple model extended to emission and absorption, showing what a difference an emitting atmosphere makes. Also very easy to see that the “IPCC logarithmic graph” is not at odds with the Beer-Lambert law.
Part Four – the effect of changing lapse rates (atmospheric temperature profile) and of overlapping the pH2O and pCO2 bands. Why surface radiation is not a mirror image of top of atmosphere radiation.
Part Five – a bit of a wrap up so far as well as an explanation of how the stratospheric temperature profile can affect “saturation”
Part Six – The Equations – the equations of radiative transfer including the plane parallel assumption and it’s nothing to do with blackbodies
Part Seven – changing the shape of the pCO2 band to see how it affects “saturation” – the wings of the band pick up the slack, in a manner of speaking
Part Eight – interesting actual absorption values of CO2 in the atmosphere from Grant Petty’s book
Part Nine – calculations of CO2 transmittance vs wavelength in the atmosphere using the 300,000 absorption lines from the HITRAN database
Part Ten – spectral measurements of radiation from the surface looking up, and from 20km up looking down, in a variety of locations, along with explanations of the characteristics
Part Eleven – Heating Rates – the heating and cooling effect of different “greenhouse” gases at different heights in the atmosphere
Part Twelve – The Curve of Growth – how absorptance increases as path length (or mass of molecules in the path) increases, and how much effect is from the “far wings” of the individual CO2 lines compared with the weaker CO2 lines
And Also -
Theory and Experiment – Atmospheric Radiation – real values of total flux and spectra compared with the theory.
The HITRAN 2008 molecular spectroscopic database, by L.S. Rothman et al, Journal of Quantitative Spectroscopy & Radiative Transfer (2009)
Note 1 – The energy absorbed from a photon, E = hv, where h = Planck’s constant, and ν= frequency. The energy of the photon has to match precisely the energy required for the change in state of the molecule.
However, due to the motion of atmospheric molecules, the actual frequency, ν’ = ν(1-ν/c), where c = speed of light, ν = frequency measured by a stationary observer.
In essence, there is now a very small (but measurable) range in frequencies (with respect to a “stationary” planet) that can be absorbed to give the precise energy level required for a transition of this molecule – and so the faster the molecules are moving, the more the absorption line “spreads out”.
Molecular speeds are given by the Maxwell-Boltzmann distribution.
Note 2 –
And the code for v0.2 reproduced here for convenience:
====== HITRAN v0.2 ===================
function [vmin vmax dv cco2 T p d vt tau ] = HITRAN_0_2()
% Absorption through atmosphere of defined temperature & pressure
% uses the HITRANS database
% v0.1 – Basic calculations to get out of the starting blocks
% v0.2 – Tidy up and allow a “slab” to be at any pressure and temperature
% with line width determined by the pressure & temperature
disp([‘—- New Run —- ‘ datestr(now) ‘ —-‘]);
% =========== Basic definitions =====================
p0=1.013e5; % std pressure for the HITRANS database in Pa
T0=296; % std temperature for HITRANS database in K
cco2=360e-6; % concentration of CO2
Na=6.02e23; % Avogradro’s number
mair=29e-3; % molar mass of air
% now define the fractional abundance of each isotopologue of CO2
isoprop=[0.98420 0.01106 0.0039471 0.000734 0.00004434 0.00000825 0.0000039573 0.00000147];
ignoreiso=false; % allows easy switching between main isotopologue only or all of them
% =========== Define wavenumber ranges ==============
vmin=640; % min wavenumber to consider (overall)
vmax=650; % max wavenumber to consider (overall)
dv=0.001; % resolution of wavenumber
% because absorption lines below the lower boundary (and above the upper
% boundary) can still affect the wavenumbers under consideration we need to
% extend the boundaries of which absorption lines to consider
rvmin=0.98; % ratio below lower boundary to consider
rvmax=1.02; % ratio above upper boundary to consider
%vt=vmin:dv:vmax; % vt = wavenumbers we will calculate optical depth for
tau=zeros(1,vtmax); % set optical depth to zero as a starting point
% =========== Read in CO2 data and select vmin-vmax data ========
% selected paramaters for complete HITRAN CO2 database already stored in file: co2.mat
% select only main isotopologue between the wavenumber boundaries
iv=find(iso==1 & v>vmin*rvmin & v<vmax*rvmax); % select upper and lower bounds
% select all isotopologues and between the wavenumber boundaries
iv=find(v>vmin*rvmin & v<vmax*rvmax); % select upper and lower bounds
% select subset of CO2 spectral data for evaluate
iso=iso(iv); v=v(iv); S=S(iv); gama=gama(iv); nair=nair(iv);
% ============= Define standard atmosphere against height ================
% reuse code from rte_0_5_0 to define std atmosphere and number of levels
% SI units used unless otherwise stated
% first a “high resolution” atmosphere
% zr = height, pr = pressure, Tr = temperature, rhor = density
Ts=296; % define surface temperature at same value as HITRAN measurements
ps=1.013e5; % define surface pressure
Ttropo=215; % define tropopause temperature
% nmv=2.079e25; % nmv x rho = total number of molecules per m^3, not yet
maxzr=50e3; % height of atmosphere
numzr=5001; % number of points used to define real atmosphere
zr=linspace(0,maxzr,numzr); % height vector from sea level to maxzr
[pr Tr rhor ztropo] = define_atmos_0_3(zr,Ts,ps,Ttropo); % function to determine (or lookup) p, T & rho
% Create “coarser resolution” atmosphere – this reduces computation
% requirements for absorption & emission of radiation
% z, p,Tinit,rho; subset of values used for RTE calcs
numz=10; % number of boundaries to consider (number of layers = numz-1)
minp=2e4; % top of atmosphere to consider in pressure (Pa)
% want to divide the atmosphere into approximately equal pressure changes
dp=(pr(1)-minp)/(numz); % finds the pressure change for each height change
zi=zeros(1,numz); % zi = lookup vector to “select” heights, pressures etc
for i=1:numz % locate each value
zi(i)=find(pr<=(pr(1)-i*dp), 1); % gets the location in the vector where
% pressure is that value
% now create the vectors of coarser resolution atmosphere
% z(1) = surface; z(numz) = TOA
% T, p, rho all need to be in the midpoint between the boundaries
% T(1) is the temperature between z(1) and z(2), etc.
z=zr(zi); % height
pb=pr(zi); % pressure at boundaries
Tb=Tr(zi); % starting temperature at boundaries
rhob=rhor(zi); % density at boundaries
% now calculate density, pressure and temperature within each layer
dz(i)=z(i+1)-z(i); % precalculate thickness of each layer
Tinit(i)=(Tb(i+1)+Tb(i))/2; % temperature in midpoint of boundary
p(i)=(pb(i+1)+pb(i))/2; % pressure in midpoint of boundary
rho(i)=(rhob(i+1)+rhob(i))/2; % density in midpoint of boundary
% but in the first version we ignore this and just calculate at std temp
% and pressure (later iterate through each layer, i
% work in cm for volumes and distances (absorption coefficients are in cm)
rhop=rhor(find(pr<=p,1)); % lookup density at pressure = p
T=Tr(find(pr<=p,1)); % lookup temperature at pressure = p
na=(Na*rhop/mair)/1e6; % number of air molecules per cm^3
Nco2=na*cco2; % number of CO2 molecules per cm^3
d=100*100; % path length in cm *** just a starting point
% ======== Iterate through each absorption line ======
% for one slab of depth, d
for j=1:ivmax % each absorption line
% v(j)=line center, S(j)=line strength, gama(j)=half width
% first do a code inefficient method – calculate the profile for each
% across the entire wavenumber range
if ignoreiso==true % if we only want to consider the main isotopologue
% ======== Iterate through each wavenumber ====
for k=1:vtmax % each wavenumber
% vt(k) is the wavenumber under consideration
title([‘Transmittance due to CO2 through ‘ num2str(round(d/100))…
‘m of atmosphere @ p= ‘ num2str(p) ‘Pa, T= ‘ num2str(T) ‘K’])
xlabel(‘Wavenumber, cm^-^1′,’FontSize’, 8)
disp([‘Range = ‘ num2str(vmin) ‘ – ‘ num2str(vmax) ‘ cm^-1, dv = ‘ num2str(dv) ‘, lines = ‘ num2str(ivmax)]);
disp([‘Number of CO2 molecules per cm^3 = ‘ num2str(Nco2) ‘, path length = ‘ num2str(d/100) ‘ m’]);
disp([‘p= ‘ num2str(p) ‘ Pa, T= ‘ num2str(T) ‘ K’]);
disp([‘—- Complete End —- ‘ datestr(now) ‘ —-‘]);
======== end of HITRAN v0.2 =====================