In Part Two, we looked at the beginnings of a very simple 1-d model for examining how the atmosphere interacts with radiation from the surface.

Simplification aids understanding so the model has some fictious gases which absorb radiation – pCO2 (pretend CO2) and pH2O (pretend water vapor). We saw that as the concentrations of pCO2 were increased we stopped seeing a change in “top of atmosphere” (TOA) flux.

That is, the “greenhouse” effect became “saturated” as the pCO2 concentration was increased.

Of course, the model was flawed. The model only included absorption of radiation by the atmosphere, with no emission. There are other over-simplifications and progressively we will try and consider them.

### Emission

Once the atmosphere can emit as well as absorb radiation the results change.

The model has been updated, and for these results is now at v3.1 (see note 5).

Here is a comparison of “no emission” vs “emission”. In each case 10 runs were carried out of different pCO2 concentrations. Each graph shows:

- TOA spectral results for runs 1, 5 and 10
- Surface spectral downward radiation (DLR or “back radiation) for the last run
- Temperature profile for the last run
- Summary graph of flux vs concentration changes for all 10 runs

No emission:

*Figure 1 – Click for a larger view*

Note that the surface downward radiation is zero. This is because the atmosphere doesn’t emit radiation in this model.

With emission:

*Figure 2 – Click for a larger image*

If you compare the spectral results you can see that in the “no emission” case the “pretend water vapor” (pH2O) band of wavenumbers 1250-1500 cm^{-1} is in “saturation”, whereas in the “emission” case, it isn’t.

This is just the natural result of an atmosphere that re-emits – and of course, gases that absorb at a given wavelength also emit at that same wavelength. (See Planck, Stefan-Boltzmann, Kirchhoff and LTE).

And a comparison of the 10 runs between the two models:

*Figure 3*

As the concentration of pCO2 increases the TOA flux reduces. Of course, we can see that this effect diminishes with more pCO2.

One interesting idea is to show these results as **radiative forcing**. *(You can find a more formal definition of the term in CO2 – An Insignificant Trace Gas? Part Seven – The Boring Numbers – although here it is not calculated according to the strict definition of allowing stratospheric adjustment)*

In essence radiative forcing is the **change** in TOA flux. When **less **flux escapes this is considered a **positive** radiative forcing. The reason is this: less flux radiated from the climate system means that less energy is leaving, which means the climate will heat (all other things being equal).

Here is the same graph expressed as radiative forcing:

*Figure 4*

If you compare it with the IPCC graph in Part Two (or Part Seven of the CO2 series) you will see it has some similarities:

*Figure 5*

Note that the actual values of radiative forcing in this model are much higher – in part because we are comparing the results of a gas with made up properties and also because we are comparing the effect from ZERO concentration of the gas vs a very high concentration.

However, one point should be clear. Using the very simple Beer-Lambert law of absorption, and the very well-known (but not so simple) Planck law of emission we find that the “radiative forcing” for changing concentrations of one gas doesn’t look like the Beer-Lambert law of absorption..

### Real World Complexity

The very simple model here – with emission – will eventually reach “saturation”.

The much more complex real-world “line by line” models eventually reach “saturation” but at much higher concentrations of CO2 than we expect to see in the climate.

This model is still a very simple model designed for education of the basics – and to allow inspection of the code.

The model has no overlaps between absorbing molecules and has no effects from the weaker lines at the edges of a band.

### Reducing Emission and “The Greenhouse Effect”

As emission of radiation is reduced due to increases in absorbing gases, with all other things being equal, the planet must warm up. It must warm up until the emission of radiation again balances the absorbed radiation (note 1).

Another way to consider the effect is to think about where the radiation to space comes from in the atmosphere. As the *opacity* of the atmosphere increases the radiation to space must be from a higher altitude. See also The Earth’s Energy Budget – Part Three.

Higher altitudes are colder and so the radiation to space is a lower value. Less radiation from the climate means the climate warms.

As the climate warms, if the lapse rate (note 3) stays the same, eventually the radiation to space – from this higher altitude – will match the absorbed solar radiation. This is how increases in radiatively active gases (aka “greenhouse” gases) affect the surface temperature (see note 4).

### The Model

The equations will be covered more thoroughly in a later article in this series.

The essence of the model is captured in this diagram for one “layer”:

*Figure 6 – One layer from the model*

The atmosphere is broken up into a number of layers. In these model runs this was set to 30 layers with a minimum pressure of 10,000 Pa (about 17km). The “boundary condition” for radiation from the surface was Planck law radiation from a surface of 300K (27°C) with an emissivity of 0.98. And the “boundary condition” for radiation from TOA was zero.

This is because we are considering “longwave radiation”. *A further complication would be to consider an absorber (like CO2) which absorbs solar radiation, but this has not been done in this model. Absorbers of “shortwave” radiation trap energy in the atmosphere rather than allowing it to be absorbed at the surface. Therefore, they affect the “in planetary balance” but don’t significantly affect the total energy absorbed.*

The transmissivity at each wavenumber for each layer was calculated. Absorptivity = 1- transmissivity (note 2).

So for each layer – and each wavenumber interval – the transmitted radiation (incident radiation x transmissivity) was calculated for each wavenumber. This was done separately for up and down radiation. The emitted radiation was calculated by the Planck formula and the emissivity (= absorptivity at that wavenumber).

The net energy absorbed as a result changed the temperature, and the model iterated through many time steps to find the final result.

With the very simple pCO2 and pH2O properties the number of timesteps made almost no difference to the TOA flux. The stratospheric temperature did vary, something for further investigation.

The wavenumber interval, dv = 5cm^{-1}, was changed to see the results. Very small changes were observed as dv changed from 5cm^{-1} to 1cm^{-1}. This is not surprising as the absorbing properties of these molecules is very simple.

### Conclusion

The model is still very simple.

The changes in TOA fluxes are significantly different for the unrealistic “no emission” case vs the “emission” case. This is to be expected.

As concentration of pCO2 increases the TOA flux reduces but by progressively smaller amounts.

When we review the results as “radiative forcing” we find that even with this simple model they resemble the IPCC “logarithmic forcing result”.

There’s more to think about with real-world gases and absorption.

Hopefully this article helps people who are trying to understand the basics a little better.

*And surely there must be mistakes in the code. **Anyone who sees anything questionable, please comment. You will probably be correct.*

People trying to do these calculations in their head, or with a pocket calculator, will be wrong. Unless they are mathematical savants. Playing the odds here.. when someone says (on this subject): “I can see that …” or “It’s clear that..” without a statement of the radiative transfer equations and boundary conditions plus their solutions to the problem – I expect that they haven’t understood the problem. And mathematical savants still need to explain to rest of us how they reached their results.

*Other articles:*

*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 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. *

### Notes

**Note 1**: Simplifications aid understanding. But in the current “climate” where many assume that climate science doesn’t understand complexity it is worth explaining a little further. Words from others neatly describe the subject of “equilibrium”:

The dominant factor determining the surface temperature of a planet is the balance between the net incoming solar radiation and outgoing thermal infra-red radiation.

That is,

radiative equilibriumcan be assumed for most practical purposes, although strictly speaking it is never achieved, as time-dependent and non-radiative processes are always present..

I.M. Vardavas & Prof. F.W. Taylor, *Radiation and Climate*, Oxford University Press (2007)

**Note 2**: Scattering of solar radiation is important, but scattering of longwave (terrestrial radiation) is very small and can be neglected. Therefore, Absorption = 1 – Transmission.

**Note 3**: Lapse rate is the temperature change with altitude – typically a reduction of 6.5K per km. This is primarily governed by “adiabatic expansion” and is affected by the amount of water vapor in the atmosphere. See Venusian Mysteries and the following articles for more.

**Note 4**: The explanations here do not include the effects of feedback. Feedback is very important, but separating different effects is important for understanding.

**Note 5**: The model still has many simplifications. The Matlab code:

It is easier to read by downloading the word doc. But pasted here below for reference. The two functions called are already documented in Part Two (and have not changed).

Matlab has many great features, but they also mean that some aspects of the code might not be clear. Feel free to ask questions.

I already see that a comment doesn’t match the code. The time step, tstep=3600*3 is 3 hrs not 12 hrs, as the comment says. However, this doesn’t alter the results.

======= v 0.3.1 ================

% 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

% 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, not yet

% used

maxzr=50e3; % height of atmosphere

numzr=2001; % 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 calculating

% pressure, temperature etc

% z, p,T,rho; subset of values used for RTE calcs

numz=30; % number of layers to consider

minp=1e4; % top of atmosphere to consider in pressure (Pa)

% want to divide the atmosphere into 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

end

% now create the vectors of coarser resolution atmosphere

z=zr(zi); % height

p=pr(zi); % pressure

T=Tr(zi); % original temperature

rho=rhor(zi); % density

for i=2:numz

dz(i)=z(i)-z(i-1); % precalculate thickness of each layer

end

% ============ 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

% work in wavenumber, cm^-1

dv=5;

v=100:dv:2500; % wavenumber (=50um – 4um)

numv=length(v);

rads=ems.*planckmv(v,Ts); % surface emissive spectral power vs wavenumber 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=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>=1250 & 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

par=zeros(1,nres); % parameter we will alter

% this section has to be changed depending on the parameter being changed

% now = CO2 concentration

par=logspace(-5,-2.5,nres); % values vary from 10^-5 (10ppm) to 10^-2.5 (3200ppm)

% par=0.01; % 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=true; % === 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) 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

end

% work out the location of subplots

if nplot==1

subr=1;subc=1; % 1 row, 1 column

elseif nplot==2

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

else

subr=3;subc=3; % 3 rows, 3 columns

end

for n=1:nres % each complete run with a new parameter to try

% — the line below has to change depending on parameter chosen

mz(2,:)=par(n);

% 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

% fixed.

% n = scenario, i = layer, j = wavenumber, k = absorber

trans=zeros(numz,numv); abso=zeros(numz,numv); % pre-allocate space

for i=2:numz % 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

end

tran(i,j)=trans; % transmissivity = 0 – 1

abso(i,j)=(1-trans)*emission; % absorptivity = emissivity = 1-transmissivity

% if emission=false, absorptivity=emissivity=0

end

end

% === 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)

% has temperature changed loop.. not yet written

% not changing pressure and density with temperature changes – only

% minor variation

tstep=3600*3; % each time step in seconds = 12 hour

for h=1:100 % main iterations to achieve equilibrium

radu=zeros(numz,numv); % initialize upward intensity at each height and wavenumber

radd=zeros(numz,numv); % initialize downward intensity at each height and wavenumber

radu(1,:)=rads; % upward surface radiation vs wavenumber

radd(end,:)=0; % 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); % zero the absorbed energy before we start

for i=2:numz % each layer

for j=1:numv % each wavenumber interval

% first calculate how much of each monochromatic ray is

% transmitted to the next layer

radu(i,j)=radu(i-1,j)*tran(i,j);

% second, calculate how much is emitted at this wavelength,

% planck function at T(i) x emissivity (=absorptivity)

% this function is spectral emissive power (pi x intensity)

radu(i,j)=radu(i,j)+abso(i,j)*3.7418e-8.*v(j)^3/(exp(v(j)*1.4388/T(i))-1);

% Change in energy = dI(v) * dv (per second)

% accumulate through each wavenumber

Eabs(i)=Eabs(i)+(radu(i,j)-radu(i-1,j))*dv;

end % each wavenumber interval

end % each layer

% Downards (have to do upward, then downward)

for i=numz-1:-1:2 % 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);

% 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(i))-1);

% accumulate energy change per second

Eabs(i)=Eabs(i)+(radu(i,j)-radu(i-1,j))*dv;

end % each wavenumber interval

dT=Eabs(i)*tstep/(cp*rho(i)*dz(i)); % change in temperature = dQ/heat capacity

T(i)=T(i)+dT; % calculate new temperature

% 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 % go through each layer

if (T(i-1)-T(i))/dz(i)>lapse % too cold, convection will readjust

T(i)=T(i-1)-(dz(i)*lapse); % adjust temperature

end

end

end

end % iterations to find equilibrium temperature

flux(n)=sum(radu(end,:)*dv); % calculate the TOA flux

% === 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

xlabel(‘Wavenumber, cm^-^1′,’FontSize’,8)

ylabel(‘W/m^2.cm^-^1′,’FontSize’,8)

title(['pCO2 = ' num2str(round(par(n)*1e6)) 'ppm, Total TOA flux= ' num2str(round(flux(n))) ' W/m^2'])

grid on

end

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′)

xlabel(‘Wavenumber, cm^-^1′,’FontSize’,8)

ylabel(‘W/m^2.cm^-^1′,’FontSize’,8)

grid on

end

if plottemp==1 % plot temperature profile vs height, if requested

plotloc=nplot-(nres>1); % get subplot location

subplot(subr,subc,plotloc),plot(T,z/1000)

title(‘Temperature vs Height’)

xlabel(‘Temperature, K’,’FontSize’,8)

ylabel(‘Height, km’,’FontSize’,8)

grid on

end

if nres>1 % produce summary plot – TOA flux vs changed parameter

subplot(subr,subc,nplot),plot(par*1e6,flux)

title(‘Summary Results’,’FontWeight’,’Bold’)

ylabel(‘TOA Flux, W/m^2′,’FontSize’,8)

xlabel(‘pCO2 ppm’,’FontSize’,8) % ==== change label for different scenarios =========

grid on

end