In Part One I said:
In the next part we will consider more advanced aspects of this subject.
However, although that is where we are heading, first I want to look at a very simple model and progressively increase its complexity. It may be illuminating.
The discussions on the Radiative Transfer Equations at the very informative blog of Judith Curry’s were very interesting. Many comments showed an appreciation for the very basics with a misunderstanding of how they fit together, for example:
IPCC declares that infrared absorption is proportional to the logarithm of GHG concentration. It is not. A logarithm might be fit to the actual curve over a small region, but it is not valid for calculations much beyond that region like IPCC’s projections. The physics governing gas absorption is the Beer-Lambert Law, which IPCC never mentions nor uses. The Beer-Lambert Law provides saturation as the gas concentration increases. IPCC’s logarithmic relation never saturates, but quickly gets silly, going out of bounds as it begins its growth to infinity.
There have also been many questions asked at this blog about the various aspects of radiation and atmosphere, so I thought some examples – via a model – might be beneficial.
The main purpose of this model is to demonstrate the effect of the basic absorption and emission processes on top of atmosphere fluxes. It is not designed to work out what troubles may lie ahead for the climate.
This model is designed to demonstrate some basic ideas about what are known as the radiative transfer equations - how radiation is absorbed and emitted as it travels from the surface to the top of atmosphere (and also from the top of atmosphere to the surface).
I aim to bridge a gap between the “blackbody shell” models and the line-by-line climate models – in a way that people can understand. Ambitious aims.
The advantage of this model should be that we can try different ideas. And the code is available for inspection.
Absorption and Beer-Lambert
The equation of absorption is known as the Beer-Lambert law – you can see more about it in CO2 – An Insignificant Trace Gas? Part Three.
In essence the Beer-Lambert law demonstrates that absorption of radiation of any wavelength is dependent on the amount of the molecule in the path of the radiation & how effective it is at absorbing at that wavelength (note 1):
Figure 1 – Beer-Lambert Law
Simple Model – v0.2 – Saturation
The simple model (v0.2) starts by creating an atmosphere. First of all, atmospheric temperatures, “created”, by definition. Second, using the temperatures plus the hydrostatic equation to calculate the pressure vs height. Third, the density (not shown) is derived.
The surface temperature is 300K.
In this first version of the model (v0.2), we introduce two radiatively-active molecules (=”greenhouse” gases):
- pH2O – pretend H2O
- pCO2 – pretend CO2
They have some vague relationship with the real molecules of similar names, but the model objective is simplicity, so of course, these molecules aren’t as complex as the real-life water vapor and CO2. And their absorption characteristics are in no way intended to be the same as the real molecules.
The mixing ratio of pH2O:
Figure 3 – water vapor mixing ratio
The pCO2 mixing ratio is constant with height, but the ppm value is varied from model run to model run.
So onto our top of atmosphere results. I’ll refer to “top of atmosphere” as TOA from now on. These results are at 20km, a fairly arbitrary choice.
The results are shown against wavenumber, which is more usual in spectroscopy, but less easy to visualize than wavelength. Wavenumber = 10,000/wavelength (when wavelength is in μm and wavenumber is in cm-1.
- 15 μm = 667 cm-1
- 10 μm = 1,000 cm-1
- 4 μm = 2,500 cm-1
The “notch” that starts to appear around 600-700 cm-1 is from pCO2 as we increase its concentration. The effect of pH2O is held constant. You can see the effect of pH2O as the “notch” around wavenumbers 1250-1500.
Note that the final graph summarizes the results – TOA Flux vs pCO2 concentration:
What we see is that eventually we reach a point where more pCO2 has no more effect – saturation!
We see that as the concentration keeps on increasing we get a reduction in the TOA flux down to a “saturated” value.
And when we plot the Summary Results we get something that looks very like the Beer-Lambert law of absorption! Nothing like the IPCC result shown in CO2 – An Insignificant Trace Gas? Part Seven – The Boring Numbers:
So, case closed then?
Flaws in the Model
This model is very simplistic. Most importantly, it only includes absorption of radiation and does not include emission.
- no emission from any layers in the atmosphere
- no dependence of line width on pressure and temperature
- no weaker “wings” of an absorption band
- no overlapping of absorption (absorption by 2 molecules in the same band)
In fact, the model only represents the Beer Lambert law of absorption and therefore is incomplete. If the atmosphere is absorbing energy from radiation, where is this energy going?
The next step is to add emission.
Another important point is that the total number of people involved in creating, testing and checking this model = 1.
Surely there will be mistakes.. I welcome comments from sharp-eyed observers pointing out the mistakes.
Part One – a bit of a re-introduction to the subject
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
And Also -
Theory and Experiment – Atmospheric Radiation – real values of total flux and spectra compared with the theory.
Note 1 – The Beer-Lambert law is written in a number of slightly different forms. The transmission of incident monochromatic radiation is reduced like this:
I(λ) = I0(λ).exp(-σ(λ)nz)
where I(λ) = transmitted radiation of wavelength λ, I0(λ) = the incident radiation, σ(λ) = “capture cross section” at wavelength λ, n = number of absorbers per unit volume, z = length of path. And this equation relies on the number of absorbers remaining constant per unit volume. The equation can be written differently to deal with the case when the density of absorbers through the path changes.
Note 2 – If there was no radiatively-active atmosphere, then the top of atmosphere flux would be the same as the surface. With an emissivity of 0.98 and a temperature of 300K, the emission of radiation from the surface = 450 W/m².
It is actually calculated in the program by numerically integrating the spectrum of surface radiation, with a result of 447 W/m² (due to not integrating the waveform with enough precision at high wavelengths/low wavenumbers).
Note 3 – The Matlab program used to create this data:
and two functions called from RTE:
==== RTE v0.2 =========
% RTE = Radiative transfer equations in atmosphere
% Objective – allow progressively more complex applications
% to help people see what happens in practice in the atmosphere
% v0.2 allow iterations of one (or more) parameter to find the TOA flux vs
% changed parameter
% Define standard atmosphere against height
% zr = height, pr = pressure, Tr = temperature, rhor = density, all in SI units
Ts=300; % define surface temperature
ps=1.013e5; % define surface pressure
nmv=2.079e25; % nmv x rho = total number of molecules per m^3
maxzr=50e3; % height of atmosphere
numzr=1000; % 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(zr,Ts,ps); % function to determine (or lookup) p, T & rho
% Consider atmosphere in a coarser resolution than used for calculation
% z, p,T,rho; subset of values used for RTE calcs
numz=20; % number of layers to consider
maxz=20e3; % height of atmosphere to consider in meters
zi=round(logspace(0.3,log10(numzr),numz)); % zi = lookup vector to “select” heights, pressures etc
% approximately logarithmically as pressure is logarithmic
% various values
lapse=6.5; % environmental lapse rate ** note potential conflict with temp profile already determined
% could be max lapse rate for convective adjustment
ems=0.98; % emissivity of surface
% work in wavenumber, cm^-1
v=100:dv:2500; % wavenumber (=50um – 4um)
rads=ems.*planckmv(v,Ts); % surface spectral intensity against frequency v
% ====== Introducing the molecules ============
% need % mixing in the atmosphere vs height, % capture cross section per
% number per frequency, pressure & temperature broadening function
nummol=2; % number of radiatively-active gases
mz=ones(nummol,numz); % initialize mixing ratios of the gases
% specific concentrations
% pH2O = pretend H2O
emax=17e-3; % max mixing ratio (surface) of 17g/kg
mz(1,:)=(ztropo-z).*emax./ztropo; % straight line reduction from surface to tropopause
mz(1,(mz(1,:)<0))=5e-6; % replace negative values with 5ppm, ie, for heights above tropopause
% pCO2 = pretend CO2
mz(2,:)=3.6e-4; % 360ppm
% absorption coefficients
k1=.3; % arbitrary pick – in m2/kg while we use rho
k2=.3; % likewise
a=zeros(nummol,length(v)); % initialize absorption coefficients
a(1,(v>=1250 & v<=1500))=k1; % wavelength dependent absorption
a(2,(v>=600 & v<=800))=k2;
% === Loop to change key parameter for which we want to see the effect
nres=10; % number of results to calculate
flux=zeros(1,nres); % TOA flux for each change in parameter
par=zeros(1,nres); % parameter we will alter
plotix=[1 3 5 8 10]; % graphs to plot
nplot=length(plotix)+1; % plot the “plotix” graphs plus the summary
% work out the location of subplots
subr=1;subc=1; % 1 row, 1 column
subr=2;subc=1; % 2 rows, 1 columns
elseif nplot==3 || nplot==4
subr=2;subc=2; % 2 rows, 2 columns
elseif nplot==5 || nplot==6
subr=3;subc=2; % 3 rows, 2 columns
subr=3;subc=3; % 3 rows, 3 columns
% parameter values to try – this section has to be changed depending on the
% CO2 concentration
for n=1:nres % each complete run with a new parameter to try
% — this line below has to change depending on parameter chosen
% === Main loops to calculate TOA spectrum & flux =====
% first we work through each layer
% then through each wavenumber
% currently calculating surface radiation absorption up to TOA
rad(1,:)=rads; % surface radiation vs wavenumber
for i=2:numz % each layer
for j=1:numv % each wavenumber interval
abs=1; % initialize the amount of absorption within the wavenumber interval
for k=1:nummol % each absorbing molecule
% for each absorber: exp(-density x mixing ratio x
% absorption coefficient x thickness of layer)
abs=abs*exp(-rho(i)*mz(k,i)*a(k,j)*(z(i)-z(i-1))); % calculate absorption
flux(n)=sum(rad(end,:)*dv); % calculate the TOA flux
% decide if and where to plot
ploc=find(plotix==n); % is this one of the results we want to plot?
if not(isempty(ploc)) % then plot. “Ploc” is the location within all the plots
subplot(subr,subc,ploc),plot(v,rad(end,:)) % plot wavenumber against TOA emissive power
title(['pCO2 = ' num2str(round(par(n)*1e6)) 'ppm, Total TOA flux= ' num2str(round(flux(n))) ' W/m^2'])
% final plot – TOA flux vs changed parameter
ylabel(‘TOA Flux, W/m^2′)
xlabel(‘pCO2 ppm’) % this x label is up for changing
====== end of RTE v0.2 =========
====== define_atmos ============
function [pr Tr rhor ztropo] = define_atmos(zr, Ts, ps)
% Calculate pressure, pr; temperature, Tr; density, rhor; from height, zr;
% surface temperature, Ts; and surface pressure, ps
% By assuming a standard temperature profile to a fixed tropopause temp
% zr is a linear vector and starts at the surface
pr=zeros(1,length(zr)); % allocate space
Tr=zeros(1,length(zr)); % allocate space
rhor=zeros(1,length(zr)); % allocate space
% first calculate temperature profile
lapset=6.5/1000; % lapse rate in K/m
Ttropo=215; % temperature of tropopause
htropo=10000; % height of tropopause
Tstrato=270; % temperature of stratopause
zstrato=50000; % height of stratopause
ztropo=((Ts-Ttropo)/lapset); % height of bottom of tropopause in m
Tr(i)=Ts-(zr(i)*lapset); % the troposphere
elseif (zr(i)>=ztropo && zr(i)<=(ztropo+htropo))
Tr(i)=Ttropo; % the tropopause
elseif (zr(i)>(ztropo+htropo) && zr(i)<zstrato)
Tr(i)=Ttropo+(zr(i)-(ztropo+htropo))*(Tstrato-Ttropo)/(zstrato-(ztropo+htropo)); % stratosphere
Tr(i)=Tstrato; % haven’t yet defined temperature above stratopause
% but this prevents an error condition
% pressure, p = ps * exp(-mg/R*integral(1/T)dz) from 0-z
mr=28.57e-3; % molar mass of air
R=8.31; % gas constant
g=9.8; % gravity constant
intzt=0; % sum the integral of dz/T in each iteration
dzr=zr(2)-zr(1); % dzr is linear so delta z is constant
pr(1)=ps; % surface condition
% density, rhor = mr.p/RT
======= end of define_atmos =================
======= planckmv ==================
function ri = planckmv( v,t )
% planck function for matrix of wavenumber, v (1/cm) and temp, t (K)
% result in spectral emissive power (pi * intensity)
% in W/(m^2.cm^-1)
[vv tt]= meshgrid(v,t);
ri = 3.7418e-8.*v.^3./(exp(vv.*1.4388./tt)-1);
============ end of planckmv =================