In Part Four we took a first look at overlapping gases. pH2O’s absorption band was changed to overlap pCO2′s absorption band. And remember that pH2O has a much higher concentration in the lower atmosphere.
For those who haven’t followed the series so far, these are fictional molecules with only a passing resemblance to the real molecules H2O and CO2. The massive complexity of real spectral absorption and emission makes it difficult for people to appreciate the key points of radiative transfer in the atmosphere.
And of course, many people don’t want to “just accept” the results of a hugely complex computer model..
The simple model results revealed some interesting points:
- With overlapping bands, increases in pCO2 still led to a reduction in TOA flux.
- With increasing pCO2, DLR (back radiation) remains constant and yet TOA flux reduces.
It’s important to understand these results, because it’s very common to see an implicit belief that the TOA results are some kind of “mirror image” of the surface downward results. They aren’t even though of course they are related.
For the results shown in Figure 7 of Part Four, here is the last TOA spectrum and below, the corresponding DLR spectrum:
The balance of energy at TOA is what determines whether the planet warms or cools. Therefore, the spectral values at the surface are not the most important for determining which gases make the most contribution to the inappropriately-named “greenhouse” effect.
And the total value of back radiation at the surface is not what determines the long term surface temperature – because it is possible to reduce the TOA flux without increasing the surface downward flux. (Note 1)
Hopefully, this simple model demonstrates those points clearly.
Just for reference I have added this model version, v0.4.0 to the notes.
Stratospheric Temperatures and “Saturation”
The model results shown in Figure 7 of Part Four show that the TOA flux continues to reduce as the pCO2 concentration increases.
There is an important point here for the ever popular theme of “saturation”.
Let’s take a look at that model again, this time up to very high concentrations of pCO2:
Figure 2 – Click for a larger image
Notice that even as the pCO2 concentration has reached 50,000ppm the TOA flux is still reducing for increasing pCO2.
Also notice the temperature profile (5th graph in figure 2) – it’s important.
Now here is a similar model run with a slightly different constraint:
Figure 3 – Click for a larger image
These results show that “saturation” is reached much sooner. Notice the temperature profile.
The intensity of radiation is dependent on the temperature of the atmosphere from where the radiation takes place.
So if we have an atmosphere that keeps reducing in temperature as we go higher, then no matter how much the concentration of a “greenhouse” gas increases, the ever-higher radiation will be from a colder temperature – and therefore, will keep reducing in intensity.
Of course, eventually the atmosphere thins out to the point where even this effect disappears.
But hopefully the basic physics behind that idea is clear.
This is why in Figure 3 where the stratospheric temperature is held constant and isothermal (all at the same temperature), the changes in TOA flux level off much sooner. No matter where in the stratosphere the atmosphere radiates from it will be at the same temperature. (See the section “Why the Lapse Rate Matters” in Part Four which is covering a very similar point).
Here is the comparison, of 20 different pCO2 concentrations, where the stratosphere was held at 215K (isothermal) and where the stratospheric temperature was allowed to change according to the radiative heating/cooling:
The temperature profile of the atmosphere does affect the “saturation” or not question by “greenhouse” gases.
Note that this model is still very simplistic – both of the gases have a fixed absorption within a band and zero outside. Real gases are much more complex and these complexities are very significant in the “saturation” question.
This article is more of a summary and consolidation so far, than any new ideas.
The next article, before covering line width issues, will cover some of the basic maths (and an explanation of the maths) behind how radiation moves through the atmosphere. At least that’s the intent at the moment.
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 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: Reducing the TOA flux = less heat leaves the planet = the planet warms; all other things being equal. More about this idea in The Earth’s Energy Budget – Part Three.
In an immediate sense the back radiation is one of the mechanisms by which the surface is at the temperature it is.
Think of the TOA flux as determining the long term temperature of the surface, and the back radiation as determining the current temperature of the surface.
And for the many who think that this means I am saying convection is unimportant, no I am not. I am explaining one effect on the surface temperature. The essence of understanding a complex subject is to be able to understand the separate effects, and then how they fit together.
Note 2: The Matlab code, v0.4.0:
The code is easiest seen by downloading the word doc, but here it is for reference:
======= v0.4.0 ======================
% 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
% v0.3 add emissivity = absorptivity ; as a function of wavelength. Also
% means that downward and upward radiation must be solved, plus iterations
% to allow temperature to change to find stable solution. Use convective
% adjustment to the lapse rate
% v0.3.1 changes the method of defining the atmosphere layers for radiation
% calculations, to have roughly constant mass for each layer
% v0.3.2 tries changing lapse rates and tropopause heights
% v0.3.3 revises element boundaries as various problems found in testing of
% v0.4.0 – introducing overlap of absorption bands
clear % empty all the variables, so previous runs can have no effect
disp(['---- New Run ---- ' datestr(now) ' ----']);
% SI units used unless otherwise stated
% ============= Define standard atmosphere against height ================
% first a “high resolution” atmosphere
% zr = height, pr = pressure, Tr = temperature, rhor = density
Ts=300; % define surface temperature
ps=1.013e5; % define surface pressure
% 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_2(zr,Ts,ps); % 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=30; % number of boundaries to consider (number of layers = numz-1)
minp=3e3; % 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
% ============ Set various values =========================
lapse=6.5e-3; % environmental lapse rate in K/m ** note potential conflict with temp profile already determined
% currently = max lapse rate for convective adjustment, but not used to
% define initial temperature profile
ems=0.98; % emissivity of surface
cp=1000; % specific heat capacity of atmosphere, J/K.kg
convadj=true; % === SET TO true === for convective adjustment to lapse rate = lapse
emission=true; % ==== SET TO true ==== for the atmosphere to emit radiation
tstep=3600*12; % fixed timestep of 1hr
nt=1000; % number of timesteps
% work in wavenumber, cm^-1
v=100:dv:2500; % wavenumber (=50um – 4um)
rads=ems.*planckmv(v,Ts); % surface emissive spectral power vs wavenumber, v
disp(['Tstep= ' num2str(tstep/3600) ' hrs , No of steps= ' num2str(nt) ', numz= ' ...
num2str(numz) ', minp= ' num2str(minp) ' Pa, Lapse= ' num2str(lapse*1e3) ' K/km']);
% ============== 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-1); % 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(1:numz-1)).*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,:)=3162e-6; % a fixed mixing ratio for pCO2
% absorption coefficients
k1=0.3; % arbitrary pick – in m2/kg while we use rho
k2=0.3; % likewise
a=zeros(nummol,length(v)); % initialize absorption coefficients
a(1,(v>=500 & v<=1500))=k1; % wavelength dependent absorption
a(2,(v>=600 & v<=800))=k2; % ” “
% ========== Scenario 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
fluxd=zeros(1,nres); % DLR for each change in parameter, not really used yet
par=zeros(1,nres); % parameter we will alter
% this section has to be changed depending on the parameter being changed
% now = pCO2 conc.
par=logspace(-5,-2.5,nres); % values vary from 10^-5 (10ppm) to 10^-2.5 (3200ppm)
% par=1; % kept for when only one value needed
% ================== Define plots required =======================
% last plot = summary but only if nres>1, ie if more than one scenario
% plot before (or last) = temperature profile, if plottemp=true
% plot before then = surface downward radiation
plottemp=false; % === SET TO true === if plot temperature profile at the end
plotdown=true; % ====SET TO true ==== if downward surface radiation required
if nres==1 % if only one scenario
plotix=1; % only one scenario graph to plot
nplot=plottemp+plotdown+1; % number of plots depends on what options chosen
else % if more than one scenario, user needs to put values below for graphs to plot
plotix=[1 round(nres/2) 8 nres]; % graphs to plot – “user” selectable
nplot=length(plotix)+plottemp+plotdown+1; % plot the “plotix” graphs plus the summary
% plus the temperature profile plus downward radiation, if required
% work out the location of subplots
subr=1;subc=1; % 1 row, 1 column
subr=1;subc=2; % 1 row, 2 columns
elseif nplot==3 || nplot==4
subr=2;subc=2; % 2 rows, 2 columns
elseif nplot==5 || nplot==6
subr=2;subc=3; % 2 rows, 3 columns
subr=3;subc=3; % 3 rows, 3 columns
for n=1:nres % each complete run with a new parameter to try
% — the line below has to change depending on parameter chosen
% to find what the stability problem is we need to store all of the
% values of T, to check the maths when it goes unstable
mz(2,:)=par(n); % this is for CO2 changes
% lapse=par(n); % this is for lapse rate changes each run
disp(['Run = ' num2str(n)]);
T=zeros(nt,numz-1); % define array to store T for each level and time step
T(1,:)=Tinit; % load temperature profile for start of scenario
% remove??? T(:,1)=repmat(Ts,nt,1); % set surface temperature as constant for each time step
% First pre-calculate the transmissivity and absorptivity of each layer
% for each wavenumber. This doesn’t change now that depth of each
% layer, number of each absorber and absorption characteristics are
% n = scenario, i = layer, j = wavenumber, k = absorber
trans=zeros(numz-1,numv); abso=zeros(numz-1,numv); % pre-allocate space
for i=1:numz-1 % each layer
for j=1:numv % each wavenumber interval
trans=1; % initialize the amount of transmission 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)
trans=trans*exp(-rho(i)*mz(k,i)*a(k,j)*dz(i)); % calculate transmission, = 1- absorption
tran(i,j)=trans; % transmissivity = 0 – 1
abso(i,j)=(1-trans)*emission; % absorptivity = emissivity = 1-transmissivity
% if emission=false, absorptivity=emissivity=0
% === Main loops to calculate TOA spectrum & flux =====
% now (v3) considering emission as well, have to find temperature stability
% first, we cycle around to confirm equilibrium temperature is reached
% second, we work through each layer
% third, through each wavenumber
% fourth, through each absorbing molecule
% currently calculating surface radiation absorption up to TOA AND
% downward radiation from TOA (at TOA = 0)
for h=2:nt % main iterations to achieve equilibrium
radu=zeros(numz,numv); % initialize upward intensity at each boundary and wavenumber
radd=zeros(numz,numv); % initialize downward intensity at each boundary and wavenumber
radu(1,:)=rads; % upward surface radiation vs wavenumber
radd(end,:)=zeros(1,numv); % downward radiation at TOA vs wavenumber
% units of radu, radd are W/m^2.cm^-1, i.e., flux per wavenumber
% h = timestep, i = layer, j = wavenumber
% Upward (have to do upward, then downward)
Eabs=zeros(numz-1); % zero the absorbed energy before we start
for i=1:numz-1 % each layer
for j=1:numv % each wavenumber interval
% first calculate how much of each monochromatic ray is
% transmitted to the next layer
% second, add emission at this wavelength:
% planck function at T(i) x emissivity (=absorptivity)
% this function is spectral emissive power (pi x intensity)
% Change in energy = dI(v) * dv (per second)
% accumulate through each wavenumber
% if the upwards radiation entering the layer is more than
% the upwards radiation leaving the layer, then a heating
end % each wavenumber interval
end % each layer
% Downwards (have to do upward, then downward)
for i=numz-1:-1:1 % each layer from the top down
for j=1:numv % each wavenumber interval
% first, calculate how much of each monochromatic ray is
% transmitted to the next layer, note that the TOA value
% is set to zero at the start
radd(i,j)=radd(i+1,j)*tran(i,j); % attentuation..
% second, calculate how much is emitted at this wavelength,
radd(i,j)=radd(i,j)+abso(i,j)*3.7418e-8.*v(j)^3/(exp(v(j)*1.4388/T(h-1,i))-1); % addition..
% accumulate energy change per second
end % each wavenumber interval
dT=Eabs(i)*tstep/(cp*rho(i)*dz(i)); % change in temperature = dQ/heat capacity
T(h,i)=T(h-1,i)+dT; % calculate new temperature
if T(h,i)>500 % Finite Element analysis problem
disp(['Terminated at n= ' num2str(n) ', h= ' num2str(h) ', z(i)= ' num2str(z(i)) ', i = ' num2str(i)]);
disp(['time = ' num2str(h*tstep/3600) ' hrs; = ' num2str(h*tstep/3600/24) ' days']);
% need a step to see how close to an equilibrium we are getting
% not yet implemented
end % each layer
% now need convective adjustment
if convadj==true % if convective adjustment chosen..
for i=2:numz-1 % go through each layer
if (T(h,i-1)-T(h,i))/dz(i)>lapse % too cold, convection will readjust
T(h,i)=T(h,i-1)-(dz(i)*lapse); % adjust temperature
end % iterations to find equilibrium temperature
flux(n)=sum(radu(end,:))*dv; % calculate the TOA flux
fluxd(n)=sum(radd(1,:))*dv; % calculate the DLR total
% === Plotting specific results =======
% 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,radu(end,:)) % plot wavenumber against TOA emissive power
title(['pCO2 @ ' num2str(round(par(n)*1e6)) 'ppm, TOA flux= ' num2str(round(flux(n)))...
' W/m^2, DLR= ' num2str(round(fluxd(n)))])
%title(['Lapse Rate ' num2str(par(n)*1000) ' K/km, Total TOA flux= ' num2str(round(flux(n))) ' W/m^2'])
end % end of each run with changed parameter to see TOA effect
if plotdown==1 % plot downward surface radiation, if requested
plotloc=nplot-plottemp-(nres>1); % get subplot location
subplot(subr,subc,plotloc),plot(v,radd(2,:)) % plot wavenumber against downward emissive power
title(['Surface Downward, W/m^2.cm^-^1, Total DLR flux= ' num2str(round(fluxd(n))) ' W/m^2'])
if plottemp==1 % plot temperature profile vs height, if requested
plotloc=nplot-(nres>1); % get subplot location
title(‘Temperature vs Height (last scenario)’)
if nres>1 % produce summary plot – TOA flux vs changed parameter
ylabel(‘TOA Flux, W/m^2′,’FontSize’,8)
xlabel(‘pCO2 concentration, ppm’,'FontSize’,8) % ==== change label for different scenarios =========
disp(['---- Complete End ---- ' datestr(now) ' ----']);