The Science of Doom

In the model simulations up until now the pCO2 band has been a constant – a fixed absorption band between the wavenumbers 600-800 cm^-1. In this article we will see what happens when the band shape changes, but the “area under the curve” stays the same.

If you are new to the series, take a look at Part Two (and probably a good idea to take a look the subsequent parts to see how the model and the results develop). pH2O and pCO2 are not the real molecules, but have a passing resemblance to them. This allowing us to see changes in top of atmosphere (TOA) flux and DLR (back radiation) as they take on different characteristics, and as the concentration of pCO2 is changed.

In Part Four we saw the effects of the pH2O band overlapping the pCO2 band. The models results in this article keep the same concept, but have widened the pH2O absorption slightly (extended to a lower wavenumber = higher wavelength).

The absorption characteristics for the four model runs:

Figure 1 – Blue is pH2O absorption, Green is pCO2 absorption

What happens is that the total band absorptance (the area under the curve) is the same for each run. The pCO2 absorption profile in Run 1 is the same as in all the previous articles in the series. Then in subsequent runs the edges of the band are widened, the “top” is made thinner, and the peak value is adjusted to keep the “total band absorptance” the same.

The code in the notes is specifically for Run 4. The code has to slightly change for each run, as I am changing the shape of the band in each run.

Each model run simulates 20 different pCO2 concentrations. The TOA spectrum for a few results is shown, along with temperature profile, absorption characteristics and the summary of TOA flux vs pCO2 concentration.

The stratospheric temperature is fixed at 215K – to understand why I’ve chosen to fix it, see the discussion in Part Five.

Run 1:

Figure 2 – Click on the image for a larger view

Run 2:

Figure 3 – Click on the image for a larger view

Run 3:

Figure 4 – Click on the image for a larger view

Run 4:

Figure 5 – Click on the image for a larger view

Now let’s compare the summary results from each of the four models:

Figure 6 – Changes in TOA flux

And also a reminder that as the TOA flux reduces, the radiative forcing increases – because a reduction in TOA flux means less energy radiated from the planet, and therefore (all other things being equal) a heating effect.

So if we plot the results as radiative forcing instead:

Figure 7 – Radiative Forcing

Effect of Band Shape

The reason for adding this complexity to the model is because real CO2 has an absorption band which is not “rectangular”. The band has a peak absorption around 667 cm^-1 (15 μm) which then falls off on either side.

In the results up until now we have seen that different effects contribute to the eventual “saturation” of increasing pCO2 . The particular conditions determine where that saturation takes place.

What these results show is that with a “rectangular” absorption profile the “saturation” takes place much sooner than with a wider band.

Perhaps that is not very surprising to many people.

If it is surprising, then for a conceptual model think about a very strong, very narrow band reaching saturation quite quickly. Whereas if you spread the same absorption over a wider band then each wavelength/wavenumber can contribute more as each wavelength takes longer to reach saturation.

I find it hard to write a good conceptual explanation. Perhaps better to suggest that readers who find the results surprising to reread Part Two to Part Five and review the results from each of the models.

Logarithmic Slopes

I also played around with a logarithmic fall-off instead of a straight line slope. The results were very similar but not identical. I applied the same normalization to the area under the curve as in the previous models.

For example, Run 8, which is similar to Run 4 (similar in that the wavenumbers for the zero points and the “flat top” are identical):

Figure 8 – Click on the image for a larger view

The value of TOA flux at the highest pCO2 concentration of Run 8 was 196.3 W/m², compared with 194.3 W/m² in Run 4.

As less absorption has been “moved out” to the “wings” of the band we would expect that the logarithmic slopes have slightly less effect than the straight line slopes.

Reducing pH2O Effectiveness

I also tried halving the absorption characteristics of pH2O to see the effect. Runs 1-4 were repeated as Runs 9-12 with the only change being the halving of the absorption coefficient of pH2O:

Figure 9

Of course, the total TOA flux is higher in all cases – as pH2O’s effectiveness has been reduced. If you compare the results with Figure 6 you see that curves are quite similar.

Conclusion

Most of the ideas that people have about “saturation” of CO2 are based on a very limited appreciation of how radiative transfer takes place.

The fact that CO2 has very strong absorption in the central part of its band has led to a false conclusion that, therefore, more CO2 can have no further effect.

As we have seen in Part Three, just considering the basic fact of emission totally changes that perspective.

Then we see the effects of emitting from different heights in the atmosphere, and how changing temperature profiles (Part Four) affects that emission.

Another area of confusion is with the idea that because water vapor has absorption across the CO2 band and because water vapor has a much higher concentration than CO2 in the lower troposphere, therefore there is nothing more for CO2 to affect. We saw in Part Four that even making the absorption coefficient of pH2O the same as pCO2  – so that the DLR (back radiation) results were swamped from pH2O – even then the TOA flux was considerably affected by increasing concentrations of pCO2. And in the real world, the real water vapor molecule has considerably less absorption than CO2 around the 15 μm / 667 cm^-1 band.

Finally we have seen how a band with a different shape causes quite different results – having a “spreading band” more like the real-world molecules also increases the effect of pCO2 at higher concentrations (less “saturation”).

Well, these are all “model results”

Which means they are the solution to some mathematical equations. If the equations are wrong, then the model is wrong (and if the model writer and his huge team made a mistake in implementation, this model is also wrong).

Yet most of the criticism of “the standard results” in the literature (and as reported by the IPCC) only shows that the critics don’t even know what the equations are.

The usual “explanations” of why the standard results are wrong fall into two categories:

1. The writers don’t state and solve the radiative transfer equations (see Part Six), but solve something else entirely. They don’t critique the standard equations, a necessary step to demonstrate the standard theory wrong. Going out on a limb, I would say that they don’t know what the standard theory is.
2. The writers create poetic words with no equations, or happily state one particular fact or figure as a breathless QED. For example, the relative concentration of water vapor vs CO2, the fact that 95% of surface radiation is absorbed within y meters by CO2 at wavelength x.

I hope that most people who have read this series have :

• a better understanding of how radiative transfer works
• a vague idea at least of what equations we would expect to see in a successful “critique”
• an appreciation of how complex the real world is and why a 2 minute calculator result is unlikely to provide insight

There may be more to come from this model. I would like to be able to include some real data on CO2 and play around with it. That might be too much work. You never know until you start.

And of course, requests and suggestions of model results you would like to see are welcome, even if no guarantees are made that the requests will be fulfilled.

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

Notes

Note 1: Sharp-eyed observers will eventually point out that I didn’t make use of the diffusivity approximation from Part Six in my model.

Note 2: The Matlab code, v0.5.0

RTE v0.5.0

The code is easiest seen by downloading the word doc, but here it is for reference:

======= v0.5.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.3.2

% v0.4.0 – introducing overlap of absorption bands

% v0.4.1 – warmer stratosphere to show “saturation” of simple bands

% v0.5.0 – introducing some band “shapes” – instead of just “rectangular”

% absorption bands. pH2O band widened to a min of 400cm^-1

clear  % empty all the variables, so previous runs can have no effect

disp(‘  ‘);

disp([‘—- New Run —- ‘ datestr(now) ‘ —-‘]);

disp(‘  ‘);

% 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

Ttropo=215; % define tropopause temperature

% nmv=2.079e25; % nmv x rho = total number of molecules per m^3, not yet

% used

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=30; % number of boundaries to consider (number of layers = numz-1)

minp=1e3; % 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

end

% 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

for i=1:numz-1

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

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

isothermal_strato=true; % ======= SET to true === for fixed and isothermal stratospheric temperature

emission=true; % ==== SET TO true ==== for the atmosphere to emit radiation

tstep=3600*12; % fixed timestep of 1hr

nt=100; % number of timesteps

% 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

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>=400 & v<=1500))=k1;  % wavelength dependent absorption

% define some parameters for the pCO2 band to make it easier to adjust

% and then “normalize” this to the original total band absorptance =

% integral(k.dv) which was 0.3 x 200cm-1 = 60

% note that v1-4 need to be multiples of dv

kint=60; % original total

v2=695;v3=705; % peak “flat top” of band

v1=400;v4=1000; % edges where band absorption coefficient falls to zero

% now optical thickness will fall off as a straight line, later we can try

% transmittance = exp(-tau) falling off as a straight line

% find the location within vector v of these points

vi1=find(v>=v1,1);vi2=find(v>=v2,1);vi3=find(v>=v3,1);vi4=find(v>=v4,1);

a(2,(v>=v2 & v<=v3))=k2;    % flat top of band

for i=vi1:vi2           % rising edge of band

a(2,i)=(v(i)-v1)/(v2-v1)*k2;

end

for i=vi3:vi4        % falling edge of band

a(2,i)=(v4-v(i))/(v4-v3)*k2;

end

knint=sum(a(2,:))*dv; % this adjusts the whole band absorption to make the

% integral the same as the original

a(2,:)=a(2,:)*kint/knint;

disp([‘v1-4= ‘ num2str(v1) ‘, ‘ num2str(v2) ‘, ‘ num2str(v3) ‘, ‘ num2str(v4) ‘, total abs for pCO2 = ‘ num2str(sum(a(2,:))*dv)]);

% plot(v,a) % to show absorption characteristics if required

% return  % if want to early for review of absorption plot

% ========== Scenario loop to change key parameter =======================

% for which we want to see the effect

%

nres=20; % 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,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

plotabs=true; % SET TO true ======= if absorption characteristics required

if nres==1   % if only one scenario

plotix=1;  % only one scenario graph to plot

nplot=plottemp+plotdown+plotabs+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=[round(nres/2) round(nres*3/4) nres]; % graphs to plot – “user” selectable

% plotix=[]; % don’t want to plot any scenarios

nplot=length(plotix)+plottemp+plotdown+plotabs+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

% 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

% fixed.

% 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

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)

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

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

% second, add emission at this wavelength:

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

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

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

% 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

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

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

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

end  % each wavenumber interval

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

if isothermal_strato && z(i)>=ztropo  % if we want to keep stratosphere at fixed temperature

T(h,i)=T(h-1,i);

else

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’]);

disp(datestr(now));

return

end

end

% 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

end

end

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

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

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

title([‘pCO2 @ ‘ num2str(round(par(n)*1e6)) ‘ppm, TOA flux= ‘ num2str(round(flux(n)))…

‘ W/m^2, DLR= ‘ num2str(round(fluxd(n)))])

% —

%subplot(subr,subc,ploc),plot(T(end,:),z(2:numz)/1000)

%title([‘Lapse Rate ‘ num2str(par(n)*1000) ‘ K/km, Total TOA flux= ‘ num2str(round(flux(n))) ‘ W/m^2’])

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

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

grid on

end

end   % end of each run with changed parameter to see TOA effect

if plotdown==true  % plot downward surface radiation, if requested

plotloc=nplot-plottemp-plotabs-(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’])

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

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

grid on

end

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

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

subplot(subr,subc,plotloc),plot(T(end,:),z(2:numz)/1000)

title(‘Temperature vs Height (last scenario)’)

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

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

grid on

end

if plotabs==true % plot show absorption characteristics if required

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

subplot(subr,subc,plotloc),plot(v,a) %

title(‘Absorption Characteristics’)

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

ylabel(‘Coefficient’,’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 concentration, ppm’,’FontSize’,8) % ==== change label for different scenarios =========

grid on

end

disp([‘—- Complete End —- ‘ datestr(now) ‘ —-‘]);

Read Full Post »

Many blogs write about over-simplifications of the radiative effects in climate. Many of these blog articles review simple explanations of how it is possible for atmospheric radiative effects to increase the surface temperature – e.g. the “blackbody shell” model.

As a result many people are confused and imagine that climate science hasn’t got past “first base” with how radiation interacts with atmospheric gases.

In any field the “over-simplified analysis” is designed to help the beginner to gain conceptual understanding of the field. Not to present the complete field of scientific endeavor.

This article will try to “bridge the gap” between the over-simplified models and the very detailed theory.

Note – it isn’t possible to cover the whole subject in one blog article and a decent treatment of radiative transfer consumes many chapters of a textbook.

There will be some maths. But I will also try to provide a non-mathematical explanation of “the maths” – or “the process”.

If you find maths daunting or incomprehensible that is understandable, but there is a lot that can be learned by trying to grasp some of the basic concepts.

Monochromatic Radiation

This means we need to treat each wavelength separately. Why? Because absorption and emission is a wavelength dependent process.

For example, here is the absorption spectrum of one part of NO2:

Figure 1

So when we consider radiation “zooming” through the atmosphere we have to take it “one wavelength” at a time.

There isn’t really any such thing as monochromatic radiation or being able to take “one wavelength at a time” – but that doesn’t stop us analyzing the problem..

A Digression on “Calculus”

How does the world of science and engineering deal with continuous change?

If a force, or a ray of radiation, or a movement of the atmosphere is a continuously changing value, how do we define it? How do we deal with it?

Calculus is the answer. This branch of mathematics allows us to deal with small changes and continuous changes and provide theorems, answers and equations.

For example, if we know something about the variable distance, s, with respect to time, t, then an equation defines the relationship between velocity, v, and these other variables:

v = ds/dt

where the “d”s at the beginning means: “the rate of change of”, so the formula means – in English:

Velocity = the rate of change of Distance with respect to Time

Generally when you see something like “da/db” it means “the rate of change of variable a with respect to variable b“.

It is also common to see Δx and δx – meaning “a small change in x”. This is different from “the continuous change of x”, but the specific rationale behind when we use “dx” and “Δx” isn’t so important for this article.

The other important area of calculus is “summing” results when again there is continuous change. If someone travels at 10 km/hr for 1 hour and then at 20km /hr for 1 hr they will have traveled 10km + 20km = 30km. That’s an easy calculation. But if velocity has continuously changed with time – how do we calculate the total distance traveled?

This means, in (harder to understand) English:

Distance = the integral of Velocity with respect to Time, between the limits of time = t1 and time = t2.

The integral is like the summation of each of the tiny distances covered in each very small time period (between t1 and t2). The integral is also often referred to as “the area under the curve”.

..end of digression

Absorption of Radiation

Let’s define a monochromatic beam of radiation, I_λ, travelling through the atmosphere:

Figure 2

We have some information we can use:

The rate of absorption of the beam of radiation as it travels through the atmosphere is proportional to the amount of absorbers at that wavelength and the ability of that absorber to absorb radiation of that wavelength

This is known as the Beer-Lambert law, and is written like this (note 1):

dI_λ = -nσI_λ .ds       [1]

which means the same thing in mathematical terms, with n=number of absorbing molecules per unit volume (corrected), and σ=capture cross-section (or effectiveness at absorbing at that wavelength), and the subscript λ indicates that this equation is only true for the radiation at this wavelength

The value σ is a material property and so constant for one gas at one wavelength at one temperature and one pressure, but varies with the temperature and pressure of the gas (see comment). The value n will depend on location in the atmosphere. If we solve this equation between two arbitrary points, s1 and s2, we get:

I_λ(s₂) = I_λ(s₁). exp [ -∫σn(s).ds ]        [2]

where the integral is between the limits s₁ and s₂

What it means in English:

The intensity of radiation at wavelength λ is reduced as it travels through the atmosphere according to the total amount of the absorber along the path. “exp” is e, or 2.718, to the power of the value in the square brackets.

If the concentration of the gas doesn’t change along the path the equation becomes a simpler version:

I_λ(s₂) = I_λ(s₁). exp [ -σn.(s₂ – s₁) ]       [3]

Optical Thickness & Transmittance

These are some important properties to understand.

Optical thickness, usually written as τ, is the property inside the exponential in equation [1].

τ = ∫σn(s).ds         [4]

where the integral is between the limits s₁ and s₂

Transmittance, usually written with a weird T symbol not available in WordPress, but with “t” here, is the amount of radiation “getting through” along the path we are interested in.

t(s₁,s₂) = exp [-τ(s₁,s₂)]       [5]

also written as:

t(s₁,s₂) = e^{-τ(s₁,s₂)}

The optical thickness is “dimensionless” as is the transmittance.

So we can rewrite equation [1] as:

I_λ(s₂) = I_λ(s₁).t(s₁,s₂)       [6]

The transmittance can be a minimum of zero – although it can never actually get to zero – and a maximum of 1. So it is simply the proportion of radiation at that wavelength which emerges through the section of atmosphere in question:

Figure 3

With optical thickness, τ = 1, transmittance, t = 0.37 – which means that 63% of the radiation is absorbed along the path and 37% is transmitted.

With optical thickness, τ = 10, t=4.5×10^-5 – that is, 45 ppm will be transmitted through the path.

A note on definitions – optical thickness is usually defined as = 0 at the top of atmosphere, where z is a maximum and a maximum at the surface, where z = 0:

Figure 4

So τ increases while z decreases, and z increases while τ decreases.

Absorptance

In the absence of scattering (note 2), absorptance, a = 1 -t.

That is, whatever doesn’t get transmitted gets absorbed.

Plane Parallel Assumption

If you refer back to Figure 2, you see that the radiation is not travelling vertically upwards, but at an angle θ to the vertical.

Using simple trigonometry, ds = dz / cos θ     [7]

It’s always an advantage if we can simplify a problem and relating everything to only the vertical height through the atmosphere helps to solve the equations.

Atmospheric properties vary much more in the vertical direction than in the horizontal direction. For example, go up 10 km and the pressure drops by a factor of 5 – from 1000 mbar to 200 mbar. But travel 10 km horizontally and the pressure will have changed by less than 1 mb. Temperature typically changes 100 times faster in the vertical direction than in the horizontal direction.

And as air density is determined by pressure and absolute temperature we can make a reasonable assumption that the density at a given height, z directly above is the same as the density at the same height when we look at an angle of 45°.

Of course, by the time we are considering an angle close to 90° – i.e., horizontal – the assumption is likely to be invalid. However, the transmissivity of the atmosphere at angles very close to the horizon is extremely low anyway, as we will see.

Therefore, making the assumption of a plane parallel atmosphere is a good approximation.

Let’s review the earlier equations using a mathematical identity that reduces “equation clutter”:

μ = cos θ      [8]

And rewrite equation [5]:

t(z₁,z₂) = exp [-τ(z₁,z₂)/μ]       [9]

Notice that the equations are now rewritten in terms of the optical thickness between two vertical heights and the angle of the radiation.

It might help to see it in graphical form – and note here that the optical thickness, τ, is for the vertical direction (otherwise the graph would make no sense):

Figure 5

This simply demonstrates that as the angle increases the radiation has to travel through more atmosphere.

So suppose the optical thickness vertically through the atmosphere, τ = 1, then for:

• a vertically travelling ray the transmittance = 0.37
• for a ray at 45° the transmittance = 0.24
• for a ray at 70° the transmittance = 0.05
• for a ray at 80° the transmittance = 0.003

Emission

Using Kirchhoff’s law, absorptivity of a material = emissivity of a material for the same wavelength and direction. For diffuse surfaces – and for gases – direction does not affect these material properties, so they are only a function of wavelength. (And for all intents and purposes, absorptance is the same term as absorptivity, and transmittance is the same term as transmissivity-see comment).

Emission of radiation at any given wavelength for a blackbody (a perfect emitter) is given by Planck’s law, which is usually annotated as B_λ(T), where T = temperature.

The absorptivity of a gas, a_λ = 1-t_λ =emissivity of a gas, ε_λ.   (Corrected)

For a very small change in monochromatic radiation due to emission:

dI_λ = ε_λB_λ(T) .ds        [10a]

Therefore:

dI_λ = nσB_λ(T) .ds        [10b]

So if now combine emission and absorption, equations 1 & 10:

dI/ds = nσ.(B_λ(T) – I_λ)    [11]

If we combine this with our definition of optical thickness, from equation [4]:

dI_λ/dτ = I_λ – B_λ(T)     [12]

which is also known as Schwarzschild’s Equation – and is the fundamental description of changes in radiation as it passes through an absorbing (and non-scattering) atmosphere.

It says, in not so easy to understand English:

The change in monochromatic radiation with respect to optical density is equal to the difference between the intensity of the radiation and the Planck (blackbody) function at the atmospheric temperature

Sorry it’s not clearer in English.

In more vernacular and less precise terms:

As radiation travels through the atmosphere, the intensity increases if the Planck blackbody emission is greater than the incoming radiation and reduces if the Planck blackbody emission is less than the incoming radiation

Solving Schwarzschild’s Equation

Notice that this important equation contains the Planck term, which is for blackbody radiation (i.e., radiation from a perfect emitter), yet the atmosphere is not a perfect emitter. We definitely haven’t assumed that the atmosphere is a blackbody and yet the Planck terms appears in the equation. It’s just how the equation “pans out”..

I mention this only because so many people have come to believe that there is some “big blackbody assumption” in climate science and they might be concerned by this term. Nothing to worry about, this has not been assumed.

Let’s find a solution to the equation if we are measuring the TOA (top of atmosphere) radiation. Refer to Figure 4:

• at TOA, z=z_m and τ=0
• at the surface, z=0 and τ = τ_m

Now some maths manipulation – skip to the end people who don’t like maths..

First we note a handy party piece:

d/dτ [Ie^-τ] = e^-τ. dI/dτ – Ie^-τ [13]

Now we multiply both sides of equation [12] by e^-τ:

e^-τ.dI_λ/dτ = e^-τ.I_λ – e^-τ.B_λ(T)       [14]

Re-arranging:

e^-τ.dI_λ/dτ – e^-τ.I_λ = – e^-τ.B_λ(T)       [14a]

And substituting “handy party piece” [13] into [14a]:

d/dτ [I_λe^-τ] = – e^-τ.B_λ(T)      [15]

Now we integrate [15] between τ=0 and τ=τ_m:

I_λ(τ_m)e^-τ_m – I_λ(0) = – ∫ B_λ(T)e^-τ dτ    [16]

where the integral is between the limits of 0 and τ_m

And re-arranging we get:

..end of maths manipulation

I_λ(0) = I_λ(τ_m)e^-τ_m + ∫ B_λ(T)e^-τ dτ    [16]

Which – for those who haven’t followed the intense maths:

I_λ(0) = I_λ(τ_m)e^-τ_m + ∫ B_λ(T)e^-τ dτ [16]

The intensity at the top of atmosphere equals..

The surface radiation attenuated by the transmittance of the atmosphere, plus..

The sum of all the contributions of atmospheric radiation – each contribution attenuated by the transmittance from that location to the top of atmosphere

Don’t worry about the maths, but it is definitely worth spending some time thinking about the words in colors – to get a conceptual understanding of how atmospheric radiation “works”.

It’s Not Over Yet – Conversion from Intensity to Flux and the Diffusivity Approximation

Remember the note about the Plane Parallel Assumption ?

Getting equations into WordPress is painful and time-consuming so a quick explanation followed by the result, especially as maths fatigue will have set in among most readers, if any made it this far.

Equation [16] is for spectral intensity. That is, one direction rather than the complete hemispherical power (flux).

To calculate spectral emissive power (flux per unit wavelength) we need to integrate the equation over one hemisphere of solid angle. We re-write equation [16] in the form of equation [9] so that the optical thickness references vertical height, z and μ, which is the cosine of the angle from the vertical. Then we integrate from μ=0 (θ=0°) to μ=1 (θ=90°).

Transmittance, t(z,0) = 2 ∫ e^-τ(z,0)/μ μ.dμ

where the integral is from 0 to 1

This equation doesn’t have an “analytical” solution, meaning we can’t rewrite it in a nice equation form without the integral. But with some clever maths that I haven’t tried to follow – but have checked – we can produce an approximation which is known as the diffusivity approximation:

2 ∫ e^-τ(z,0)/μ μ.dμ ≈ e^-τ/μ’

Where μ’ is the “effective angle” which produces a close approximation to the actual answer without needing to integrate across all angles (for each wavelength). The best value of μ’ = 0.6.

Here is the calculated integral (left side of the equation) vs the approximation with μ’ = 0.6, as a function of optical thickness, τ, demonstrating the usefulness of the approximation:

Figure 6

There are some other refinements needed, for example, the reflection of atmospheric radiation for a surface emissivity < 1, which is then attenuated by the absorptance of the atmosphere before contributing the TOA measurement. But these factors can all be introduced into the equations.

Full Color Solution

What we have produced so far is a solution for each monochromatic wavelength, λ.

Also, we haven’t explicitly stated the fact that the optical thickness depends on the concentration and “capture cross section” of each absorber for that wavelength. So the optical thickness, and transmittance, for each height requires combining the effects of each active molecule.

The solution for flux, W/m², requires integrating the equations across all wavelengths.

Wow. Conceptually straightforward. But computationally very expensive – check out Figure 1 – the absorption characteristics of each radiatively-active gas vary significantly with wavelength.

Conclusion

The equation for radiative transfer is commonly known, (in differential form) as Schwarzschild’s Equation. It relies on fundamental physics.

To solve the equation requires some maths.

To solve the equation in practical terms the plane parallel assumption is used. This relies on the fact that variations in temperature and pressure (and therefore density) are negligible in the horizontal direction compared with the vertical direction.

The equation could be solved without this plane parallel assumption, but the variations horizontally in pressure and temperature are so slight that the same result would be obtained, unless extremely high quality data on temperature, pressure, density and concentration of absorbers was available.

To solve the equation in practical terms we need to know:

• the temperature (vs height) in the atmosphere
• the concentration of each absorber vs height
• the absorption characteristics of each absorber vs wavelength

In any practical field, the “proof of the pudding is in the eating”, and so take a look at Theory and Experiment – Atmospheric Radiation – where theoretical and practical results are compared.

And lastly, the Stefan-Boltzmann equation, correct and accurate though it is (check out Planck, Stefan-Boltzmann, Kirchhoff and LTE) is not used in the actual equations of radiative transfer in the atmosphere. Nor is any assumption of “unrealistic blackbodies”.

I only note these last points due to the high quantity (but not high quality), of blog articles and comments demonstrating the writers haven’t actually read a textbook on the subject, but still feel qualified to pass judgement on this field of scientific endeavor.

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 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 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: There are many formulations of the Beer-Lambert law and even much dispute about who exactly the law should be attributed to.

Other formulations include using the density of the gas and a matching coefficient for the effectiveness of the gas at absorbing.

Note 2: When considering solar radiation (shortwave), scattering is important. When considering terrestrial radiation (longwave), scattering can be neglected. In this article, we will ignore scattering, so the results will be appropriate for longwave but not correct for shortwave.

Read Full Post »

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:

1. With overlapping bands, increases in pCO2 still led to a reduction in TOA flux.
2. 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:

Figure 1

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:

Log plot

Linear plot

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.

Conclusion

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.

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

Notes

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:

RTE 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.3.2

% v0.4.0 – introducing overlap of absorption bands

clear  % empty all the variables, so previous runs can have no effect

disp(‘  ‘);

disp([‘—- New Run —- ‘ datestr(now) ‘ —-‘]);

disp(‘  ‘);

% 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

% used

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

end

% 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

for i=1:numz-1

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

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

tstep=3600*12; % fixed timestep of 1hr

nt=1000; % number of timesteps

% 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

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

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

% 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

% fixed.

% 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

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)

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

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

% second, add emission at this wavelength:

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

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

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

% 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

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

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

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

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’]);

disp(datestr(now));

return

end

% 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

end

end

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

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

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

title([‘pCO2 @ ‘ num2str(round(par(n)*1e6)) ‘ppm, TOA flux= ‘ num2str(round(flux(n)))…

‘ W/m^2, DLR= ‘ num2str(round(fluxd(n)))])

% —

%subplot(subr,subc,ploc),plot(T(end,:),z(2:numz)/1000)

%title([‘Lapse Rate ‘ num2str(par(n)*1000) ‘ K/km, Total TOA flux= ‘ num2str(round(flux(n))) ‘ W/m^2’])

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

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

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, Total DLR flux= ‘ num2str(round(fluxd(n))) ‘ W/m^2’])

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(end,:),z(2:numz)/1000)

title(‘Temperature vs Height (last scenario)’)

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 concentration, ppm’,’FontSize’,8) % ==== change label for different scenarios =========

grid on

end

disp([‘—- Complete End —- ‘ datestr(now) ‘ —-‘]);

Read Full Post »

In Part Three we saw some results of absorption and emission in the atmosphere but later I commented that the model had some flaws.

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.

Overlapping Bands

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:

Figure 5

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:

Figure 6

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.

Conclusion

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.

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

Notes

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

RTE 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

% v0.3.2

clear  % empty all the variables, so previous runs can have no effect

disp(‘  ‘);

disp([‘—- New Run —- ‘ datestr(now) ‘ —-‘]);

disp(‘  ‘);

% 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

% used

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

end

% 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

for i=1:numz-1

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

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

tstep=3600*12; % fixed timestep of 1hr

nt=1000; % number of timesteps

% 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

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

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

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

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

% fixed.

% 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

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)

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

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

% second, add emission at this wavelength:

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

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

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

% 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

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

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

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

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’]);

disp(datestr(now));

return

end

% 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

end

end

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

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

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

title([‘pCO2 conc= ‘ num2str(round(par(n)*1e6)) ‘ppm, TOA flux= ‘ num2str(round(flux(n)))…

‘ W/m^2, DLR= ‘ num2str(round(fluxd(n)))])

% —

% subplot(subr,subc,ploc),plot(T(end,:),z(2:numz)/1000)

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

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, Total DLR flux= ‘ num2str(round(fluxd(n))) ‘ W/m^2’])

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(end,:),z(2:numz)/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,flux)

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

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

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

grid on

end

disp([‘—- Complete End —- ‘ datestr(now) ‘ —-‘]);

Read Full Post »