Take a look at Part Three for more about the model.
The Model Flaws and Update
The model (v0.3.1) had poorly defined boundaries and as a result the downward flux through a layer of the atmosphere was affecting the temperature of an adjacent layer.
I noticed the problems when looking at the model stratospheric temperature profile (the upper atmosphere). According to theory, with no ozone the stratospheric temperature shouldn’t increase from the tropopause. Yet in my model it did. So I ran some stability tests:
- different number of layers
- different top of atmosphere height/pressure
- longer runs
- shorter and longer time steps
Flaws kept appearing, including instabilities.
Anyway, the model is now fixed (according to the army of Science of Doom testers). v0.3.3 is logged in the notes.
Because the flaw was in the stratosphere, fixing the flaw (luckily) had minimal effect on the TOA flux results previously reported. Here are the new results:
Figure 1 – Click on image for a larger view
Note: ppm concentration in “Summary Results” has an incorrect legend. It is not ppm, but just a mixing ratio. The highest value is just over 3×10-3 , i.e., just over 3000 ppm.
If we compare with Figure 2 in Part Three, we see the results are similar.
Here are some stability test results after the model was fixed, trying different time steps for the same total time, and comparing the final temperature profile as well as TOA flux:
Figure 2 – Click on image for a larger view
Note: there is an error in the graph title. The model time period was 10,000 hrs (60 weeks).
And trying a fixed time step with different number of timesteps (so the total time is variable):
Figure 3 – Click on image for a larger view
So with the earlier mistakes out of the way..
Why the Lapse Rate Matters
Here are the results from 8 runs with fixed pCO2 concentration (at the highest value from earlier runs) and different lapse rates (note 1). Only 4 of the individual results are shown, and the temperature profile is from the last run:
Figure 4 – Click on image for a larger view
The first graph shows that the TOA flux is equal to the surface flux – and the spectrum shows no characteristic “notches” in it. Yet there are 2 highly absorbing gases (pH2O and pCO2) present. How can this be?
The radiatively-active gases absorb and also emit. If the emission is from a location in the atmosphere which is at the same temperature as the source of the original radiation then, in simple terms, “the amount taken out (absorbed)” = “the amount put back in (emitted)”.
The lapse rate increases in each of the following 3 graphs in figure 4, meaning that the atmosphere gets colder at any given height. Therefore, the emitted radiation from a given height will be from a colder gas and therefore will be of a lower intensity.
Question – if the lapse rate in the last scenario run is 10 K/km what do we learn from the 5th graph, which shows temperature vs height?
What the above model runs show is that changes in the lapse rate affect the inappropriately-named “greenhouse” effect = the difference between the surface radiation and the TOA radiation.
For example, if there is more water vapor in the lower atmosphere the lapse rate will reduce. As a result the TOA flux will increase. And, as already explained in Part Three (under “Reducing Emission and “The Greenhouse Effect”), this will have the effect of reducing the surface temperature – all other things being equal. Of course, more water vapor will also change the atmospheric absorption so all other things aren’t equal.
A common area of concern for people trying to understand the effect of absorbing gases in the atmosphere is this:
How can CO2 have any impact when water vapor has a much higher concentration in the lower atmosphere, has such a high absorption and overlaps the CO2 band?
The model results might be interesting here. The model is still very simplistic. Here are the new absorption characteristics:
For no particular reason the absorption coefficients are the same, and the pH2O absorption now completely encompasses the pCO2 band.
As in Part Two, the concentration of pH2O in the lower atmosphere is much higher than pCO2, even for the highest concentration of pCO2 simulated:
Here are the results of the simulation runs:
Figure 7 – Click on the image for a larger view
This should be of interest. pH2O has a much higher concentration in the lower atmosphere. Yet, as the concentration of pCO2 increases the TOA flux is affected significantly. Of course, not nearly as much as when pCO2 alone affected 600-800 cm-1, and the concentration has to reach a higher value than before to affect the TOA flux.
Notice as well that the DLR, or back radiation, is constant in each of the graphs in the last figure. How can this be?
Plenty of food for thought.
Last point for reference – these are fictional molecules, created for the purpose of illustrating the effect of absorption and emission through the atmosphere. They have vaguely similar characteristics to the real molecules H2O and CO2, but are far from identical.
The simplistic model demonstrates that lapse rate plays an important part in the effect on top of atmosphere fluxes.
The model also demonstrates that even with a higher concentration absorber in the lower atmosphere, an absorber higher up in the atmosphere can have a significant effect.
These kind of results are not easy to calculate in your head. That’s just because the equations of radiative transfer, although well-known, are not linear. And not many people can calculate summations – in their head – across multiple changing variables when one of the key terms is e-x.
The model is still quite simple, not including effects like line width and its dependence on pressure and temperature. And not including the hideous complexity that is stored in the HITRANS database.
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 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: 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 2: The Matlab code, v0.3.3
The code is easiest seen by downloading the word doc, but here it is for reference:
======= v0.3.3 ======================
% 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
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=1e4; % 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*1000) ' 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,:)=3e-3; % 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>=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
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 = CO2 concentration
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=true; % === SET TO true === if plot temperature profile at the end
plotdown=false; % ====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 4 7 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
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 conc= ' num2str(round(par(n)*1e6)) 'ppm, TOA flux= ' num2str(round(flux(n)))...
' W/m^2, DLR= ' num2str(round(fluxd(n)))])
% title(['No of time steps = ' num2str(par(n)) ', Total TOA flux= ' num2str(round(flux(n))) ' W/m^2'])
% xlabel(‘Temperature, K’,’FontSize’,8)
% ylabel(‘Height, km’,’FontSize’,8)
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’)
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) ' ----']);