Feeds:
Posts
Comments

Archive for January, 2011

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

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

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

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

Emission

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

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

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

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

No emission:

Figure 1 – Click for a larger view

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

With emission:

Figure 2 – Click for a larger image

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

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

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

Figure 3

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

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

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

Here is the same graph expressed as radiative forcing:

Figure 4

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

Figure 5

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

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

Real World Complexity

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

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

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

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

Reducing Emission and “The Greenhouse Effect”

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

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

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

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

The Model

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

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

Figure 6 – One layer from the model

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

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

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

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

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

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

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

Conclusion

The model is still very simple.

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

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

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

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

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

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

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

Other articles:

Part One – a bit of a re-introduction to the subject

Part Two – introducing a simple model, with molecules pH2O and pCO2 to demonstrate some basic effects in the atmosphere. This part – absorption only

Part Four – the effect of changing lapse rates (atmospheric temperature profile) and of overlapping the pH2O and pCO2 bands. Why surface radiation is not a mirror image of top of atmosphere radiation.

Part Five – a bit of a wrap up so far as well as an explanation of how the stratospheric temperature profile can affect “saturation”

Part Six – The Equations – the equations of radiative transfer including the plane parallel assumption and it’s nothing to do with blackbodies

Part Seven – changing the shape of the pCO2 band to see how it affects “saturation” – the wings of the band pick up the slack, in a manner of speaking

And Also -

Theory and Experiment – Atmospheric Radiation – real values of total flux and spectra compared with the theory.


Notes

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

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

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

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

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

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

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

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

RTE v0.3.1

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

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

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

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

% RTE = Radiative transfer equations in atmosphere

% Objective – allow progressively more complex applications

% to help people see what happens in practice in the atmosphere

% v0.2 allow iterations of one (or more) parameter to find the TOA flux vs

% changed parameter

% v0.3 add emissivity = absorptivity ; as a function of wavelength. Also

% means that downward and upward radiation must be solved, plus iterations

% to allow temperature to change to find stable solution. Use convective

% adjustment to the lapse rate

% v0.3.1 changes the method of defining the atmosphere layers for radiation

% calculations, to have roughly constant mass for each layer

% Define standard atmosphere against height

% zr = height, pr = pressure, Tr = temperature, rhor = density, all in SI units

Ts=300; % define surface temperature

ps=1.013e5; % define surface pressure

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

% used

maxzr=50e3; % height of atmosphere

numzr=2001; % number of points used to define real atmosphere

zr=linspace(0,maxzr,numzr); % height vector from sea level to maxzr

[pr Tr rhor ztropo] = define_atmos(zr,Ts,ps); % function to determine (or lookup) p, T & rho

% Consider atmosphere in a coarser resolution than used for calculating

% pressure, temperature etc

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

numz=30; % number of layers to consider

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

% want to divide the atmosphere into equal pressure changes

dp=(pr(1)-minp)/(numz); % finds the pressure change for each height change

zi=zeros(1,numz); % zi = lookup vector to “select” heights, pressures etc

for i=1:numz % locate each value

zi(i)=find(pr<=(pr(1)-i*dp), 1); % gets the location in the vector where

% pressure is that value

end

% now create the vectors of coarser resolution atmosphere

z=zr(zi);   % height

p=pr(zi);   % pressure

T=Tr(zi);   % original temperature

rho=rhor(zi);  % density

for i=2:numz

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

end

% ============ Set various values =========================

lapse=6.5e-3; % environmental lapse rate in K/m ** note potential conflict with temp profile already determined

% currently = max lapse rate for convective adjustment, but not used to

% define initial temperature profile

ems=0.98; % emissivity of surface

cp=1000; % specific heat capacity of atmosphere, J/K.kg

convadj=true; % === SET TO true === for convective adjustment to lapse rate = lapse

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

% work in wavenumber, cm^-1

dv=5;

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

numv=length(v);

rads=ems.*planckmv(v,Ts); % surface emissive spectral power vs wavenumber v

% ============== Introducing the molecules ==============================

% need % mixing in the atmosphere vs height, % capture cross section per

% number per frequency, pressure & temperature broadening function

nummol=2; % number of radiatively-active gases

mz=ones(nummol,numz); % initialize mixing ratios of the gases

% specific concentrations

% pH2O = pretend H2O

emax=17e-3; % max mixing ratio (surface) of 17g/kg

mz(1,:)=(ztropo-z).*emax./ztropo; % straight line reduction from surface to tropopause

mz(1,(mz(1,:)<0))=5e-6; % replace negative values with 5ppm, ie, for heights above tropopause

% pCO2 = pretend CO2

mz(2,:)=3.6e-4; % 360ppm

% absorption coefficients

k1=0.3; % arbitrary pick – in m2/kg while we use rho

k2=0.3; % likewise

a=zeros(nummol,length(v)); % initialize absorption coefficients

a(1,(v>=1250 & v<=1500))=k1;  % wavelength dependent absorption

a(2,(v>=600 & v<=800))=k2;    % ” “

% === Scenario loop to change key parameter for which we want to see the effect

%

nres=10; % number of results to calculate

flux=zeros(1,nres); % TOA flux for each change in parameter

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

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

% now = CO2 concentration

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

% par=0.01; % kept for when only one value needed

% ================== Define plots required =======================

% last plot = summary but only if nres>1, ie if more than one scenario

% plot before (or last) = temperature profile, if plottemp=true

% plot before then = surface downward radiation

plottemp=true; % === SET TO true === if plot temperature profile at the end

plotdown=true; % ====SET TO true ==== if downward surface radiation required

if nres==1   % if only one scenario

plotix=1;  % only one scenario graph to plot

nplot=plottemp+plotdown+1;  % number of plots depends on what options chosen

else     % if more than one scenario, user needs to put values below for graphs to plot

plotix=[1 round(nres/2) nres]; % graphs to plot – “user” selectable

nplot=length(plotix)+plottemp+plotdown+1; % plot the “plotix” graphs plus the summary

% plus the temperature profile plus downward radiation, if required

end

% work out the location of subplots

if nplot==1

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

elseif nplot==2

subr=1;subc=2;  % 1 row, 2 columns

elseif nplot==3 || nplot==4

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

elseif nplot==5 || nplot==6

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

else

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

end

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

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

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

% First pre-calculate the transmissivity and absorptivity of each layer

% for each wavenumber. This doesn’t change now that depth of each

% layer, number of each absorber and absorption characteristics are

% fixed.

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

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

for i=2:numz   % each layer

for j=1:numv   % each wavenumber interval

trans=1;  % initialize the amount of transmission within the wavenumber interval

for k=1:nummol  % each absorbing molecule

% for each absorber: exp(-density x mixing ratio x

% absorption coefficient x thickness of layer)

trans=trans*exp(-rho(i)*mz(k,i)*a(k,j)*dz(i)); % calculate transmission, = 1- absorption

end

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

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

% if emission=false, absorptivity=emissivity=0

end

end

% === Main loops to calculate TOA spectrum & flux =====

% now (v3) considering emission as well, have to find temperature stability

% first, we cycle around to confirm equilibrium temperature is reached

% second, we work through each layer

% third, through each wavenumber

% fourth, through each absorbing molecule

% currently calculating surface radiation absorption up to TOA AND

% downward radiation from TOA (at TOA = 0)

% has temperature changed loop.. not yet written

% not changing pressure and density with temperature changes – only

% minor variation

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

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

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

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

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

radd(end,:)=0;  % downward radiation at TOA vs wavenumber

% units of radu, radd are W/m^2.cm^-1, i.e., flux per wavenumber

% h = timestep, i = layer, j = wavenumber

% Upward (have to do upward, then downward)

Eabs=zeros(numz); % zero the absorbed energy before we start

for i=2:numz   % each layer

for j=1:numv   % each wavenumber interval

% first calculate how much of each monochromatic ray is

% transmitted to the next layer

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

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

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

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

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

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

% accumulate through each wavenumber

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

end  % each wavenumber interval

end  % each layer

% Downards (have to do upward, then downward)

for i=numz-1:-1:2   % each layer from the top down

for j=1:numv   % each wavenumber interval

% first, calculate how much of each monochromatic ray is

% transmitted to the next layer, note that the TOA value

% is set to zero at the start

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

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

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

% accumulate energy change per second

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

end  % each wavenumber interval

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

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

% need a step to see how close to an equilibrium we are getting

% not yet implemented

end  % each layer

% now need convective adjustment

if convadj==true  % if convective adjustment chosen..

for i=2:numz  % go through each layer

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

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

end

end

end

end   % iterations to find equilibrium temperature

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

% === Plotting specific results =======

% Decide if and where to plot

ploc=find(plotix==n);  % is this one of the results we want to plot?

if not(isempty(ploc))  % then plot. “Ploc” is the location within all the plots

subplot(subr,subc,ploc),plot(v,radu(end,:))  % plot wavenumber against TOA emissive power

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

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

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

grid on

end

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

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

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

subplot(subr,subc,plotloc),plot(v,radd(2,:))  % plot wavenumber against downward emissive power

title(‘Surface Downward, W/m^2.cm^-^1′)

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

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

grid on

end

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

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

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

title(‘Temperature vs Height’)

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

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

grid on

end

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

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

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

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

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

grid on

end

Read Full Post »

In Part One I said:

In the next part we will consider more advanced aspects of this subject.

However, although that is where we are heading, first I want to look at a very simple model and progressively increase its complexity. It may be illuminating.

The discussions on the Radiative Transfer Equations at the very informative blog of Judith Curry’s were very interesting. Many comments showed an appreciation for the very basics with a misunderstanding of how they fit together, for example:

IPCC declares that infrared absorption is proportional to the logarithm of GHG concentration. It is not. A logarithm might be fit to the actual curve over a small region, but it is not valid for calculations much beyond that region like IPCC’s projections. The physics governing gas absorption is the Beer-Lambert Law, which IPCC never mentions nor uses. The Beer-Lambert Law provides saturation as the gas concentration increases. IPCC’s logarithmic relation never saturates, but quickly gets silly, going out of bounds as it begins its growth to infinity.

There have also been many questions asked at this blog about the various aspects of radiation and atmosphere, so I thought some examples – via a model – might be beneficial.

The Model

The main purpose of this model is to demonstrate the effect of the basic absorption and emission processes on top of atmosphere fluxes. It is not designed to work out what troubles may lie ahead for the climate.

This model is designed to demonstrate some basic ideas about what are known as the radiative transfer equations - how radiation is absorbed and emitted as it travels from the surface to the top of atmosphere (and also from the top of atmosphere to the surface).

I aim to bridge a gap between the “blackbody shell” models and the line-by-line climate models – in a way that people can understand. Ambitious aims.

The advantage of this model should be that we can try different ideas. And the code is available for inspection.

Absorption and Beer-Lambert

The equation of absorption is known as the Beer-Lambert law – you can see more about it in CO2 – An Insignificant Trace Gas? Part Three.

In essence the Beer-Lambert law demonstrates that absorption of radiation of any wavelength is dependent on the amount of the molecule in the path of the radiation & how effective it is at absorbing at that wavelength (note 1):

Figure 1 – Beer-Lambert Law

Simple Model – v0.2 – Saturation

The simple model (v0.2) starts by creating an atmosphere. First of all, atmospheric temperatures, “created”, by definition. Second, using the temperatures plus the hydrostatic equation to calculate the pressure vs height. Third, the density (not shown) is derived.

The surface temperature is 300K.

Figure 2

In this first version of the model (v0.2), we introduce two radiatively-active molecules (=”greenhouse” gases):

  • pH2O – pretend H2O
  • pCO2 – pretend CO2

They have some vague relationship with the real molecules of similar names, but the model objective is simplicity, so of course, these molecules aren’t as complex as the real-life water vapor and CO2. And their absorption characteristics are in no way intended to be the same as the real molecules.

The mixing ratio of pH2O:

Figure 3 – water vapor mixing ratio

The pCO2 mixing ratio is constant with height, but the ppm value is varied from model run to model run.

So onto our top of atmosphere results. I’ll refer to “top of atmosphere” as TOA from now on. These results are at 20km, a fairly arbitrary choice.

The results are shown against wavenumber, which is more usual in spectroscopy, but less easy to visualize than wavelength. Wavenumber = 10,000/wavelength (when wavelength is in μm and wavenumber is in cm-1.

  • 15 μm = 667 cm-1
  • 10 μm = 1,000 cm-1
  • 4 μm = 2,500 cm-1

The “notch” that starts to appear around 600-700 cm-1 is from pCO2 as we increase its concentration. The effect of pH2O is held constant. You can see the effect of pH2O as the “notch” around wavenumbers 1250-1500.

Note that the final graph summarizes the results – TOA Flux vs pCO2 concentration:

Figure 4

What we see is that eventually we reach a point where more pCO2 has no more effect – saturation!

We see that as the concentration keeps on increasing we get a reduction in the TOA flux down to a “saturated” value.

And when we plot the Summary Results we get something that looks very like the Beer-Lambert law of absorption! Nothing like the IPCC result shown in CO2 – An Insignificant Trace Gas? Part Seven – The Boring Numbers:

Radiative Forcing vs CO2 concentration, Myhre et al (1998)

Radiative Forcing vs CO2 concentration, Myhre et al (1998)

Figure 5

So, case closed then?

Flaws in the Model

This model is very simplistic. Most importantly, it only includes absorption of radiation and does not include emission.

Simplifications:

  • no emission from any layers in the atmosphere
  • no dependence of line width on pressure and temperature
  • no weaker “wings” of an absorption band
  • no overlapping of absorption (absorption by 2 molecules in the same band)

In fact, the model only represents the Beer Lambert law of absorption and therefore is incomplete. If the atmosphere is absorbing energy from radiation, where is this energy going?

The next step is to add emission.

Another important point is that the total number of people involved in creating, testing and checking this model = 1.

Surely there will be mistakes.. I welcome comments from sharp-eyed observers pointing out the mistakes.

Other articles:

Part One – a bit of a re-introduction to the subject

Part Three – the simple model extended to emission and absorption, showing what a difference an emitting atmosphere makes. Also very easy to see that the “IPCC logarithmic graph” is not at odds with the Beer-Lambert law.

Part Four – the effect of changing lapse rates (atmospheric temperature profile) and of overlapping the pH2O and pCO2 bands. Why surface radiation is not a mirror image of top of atmosphere radiation.

Part Five – a bit of a wrap up so far as well as an explanation of how the stratospheric temperature profile can affect “saturation”

Part Six – The Equations – the equations of radiative transfer including the plane parallel assumption and it’s nothing to do with blackbodies

Part Seven – changing the shape of the pCO2 band to see how it affects “saturation” – the wings of the band pick up the slack, in a manner of speaking

And Also -

Theory and Experiment – Atmospheric Radiation – real values of total flux and spectra compared with the theory.


Notes

Note 1 – The Beer-Lambert law is written in a number of slightly different forms. The transmission of incident monochromatic radiation is reduced like this:

I(λ) = I0(λ).exp(-σ(λ)nz)

where I(λ) = transmitted radiation of wavelength λ, I0(λ) = the incident radiation, σ(λ) = “capture cross section” at wavelength λ, n = number of absorbers per unit volume, z = length of path. And this equation relies on the number of absorbers remaining constant per unit volume. The equation can be written differently to deal with the case when the density of absorbers through the path changes.

Note 2 – If there was no radiatively-active atmosphere, then the top of atmosphere flux would be the same as the surface. With an emissivity of 0.98 and a temperature of 300K, the emission of radiation from the surface = 450 W/m².

It is actually calculated in the program by numerically integrating the spectrum of surface radiation, with a result of 447 W/m² (due to not integrating the waveform with enough precision at high wavelengths/low wavenumbers).

Note 3 – The Matlab program used to create this data:

RTE v0.2

and two functions called from RTE:

define_atmos

planckmv

==== RTE v0.2 =========

% RTE = Radiative transfer equations in atmosphere

% Objective – allow progressively more complex applications

% to help people see what happens in practice in the atmosphere

% v0.2 allow iterations of one (or more) parameter to find the TOA flux vs

% changed parameter

% Define standard atmosphere against height

% zr = height, pr = pressure, Tr = temperature, rhor = density, all in SI units

Ts=300; % define surface temperature

ps=1.013e5; % define surface pressure

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

maxzr=50e3; % height of atmosphere

numzr=1000; % number of points used to define real atmosphere

zr=linspace(0,maxzr,numzr); % height vector from sea level to maxzr

[pr Tr rhor ztropo] = define_atmos(zr,Ts,ps); % function to determine (or lookup) p, T & rho

% Consider atmosphere in a coarser resolution than used for calculation

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

numz=20; % number of layers to consider

maxz=20e3; % height of atmosphere to consider in meters

zi=round(logspace(0.3,log10(numzr),numz)); % zi = lookup vector to “select” heights, pressures etc

% approximately logarithmically as pressure is logarithmic

z=zr(zi);

p=pr(zi);

T=Tr(zi);

rho=rhor(zi);

% various values

lapse=6.5; % environmental lapse rate ** note potential conflict with temp profile already determined

% could be max lapse rate for convective adjustment

ems=0.98; % emissivity of surface

% work in wavenumber, cm^-1

dv=5;

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

numv=length(v);

rads=ems.*planckmv(v,Ts); % surface spectral intensity against frequency v

% ====== Introducing the molecules ============

% need % mixing in the atmosphere vs height, % capture cross section per

% number per frequency, pressure & temperature broadening function

nummol=2; % number of radiatively-active gases

mz=ones(nummol,numz); % initialize mixing ratios of the gases

% specific concentrations

% pH2O = pretend H2O

emax=17e-3; % max mixing ratio (surface) of 17g/kg

mz(1,:)=(ztropo-z).*emax./ztropo; % straight line reduction from surface to tropopause

mz(1,(mz(1,:)<0))=5e-6; % replace negative values with 5ppm, ie, for heights above tropopause

% pCO2 = pretend CO2

mz(2,:)=3.6e-4; % 360ppm

% absorption coefficients

k1=.3; % arbitrary pick – in m2/kg while we use rho

k2=.3; % likewise

a=zeros(nummol,length(v)); % initialize absorption coefficients

a(1,(v>=1250 & v<=1500))=k1;  % wavelength dependent absorption

a(2,(v>=600 & v<=800))=k2;

% === Loop to change key parameter for which we want to see the effect

%

nres=10; % number of results to calculate

flux=zeros(1,nres); % TOA flux for each change in parameter

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

plotix=[1 3 5 8 10]; % graphs to plot

nplot=length(plotix)+1; % plot the “plotix” graphs plus the summary

% work out the location of subplots

if nplot==1

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

elseif nplot==2

subr=2;subc=1;  % 2 rows, 1 columns

elseif nplot==3 || nplot==4

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

elseif nplot==5 || nplot==6

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

else

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

end

% parameter values to try – this section has to be changed depending on the

% run

% CO2 concentration

par=logspace(-5,-2.5,nres);

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

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

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

% === Main loops to calculate TOA spectrum & flux =====

% first we work through each layer

% then through each wavenumber

% currently calculating surface radiation absorption up to TOA

rad=zeros(numz,numv);

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

for i=2:numz   % each layer

for j=1:numv   % each wavenumber interval

abs=1;  % initialize the amount of absorption within the wavenumber interval

for k=1:nummol  % each absorbing molecule

% for each absorber: exp(-density x mixing ratio x

% absorption coefficient x thickness of layer)

abs=abs*exp(-rho(i)*mz(k,i)*a(k,j)*(z(i)-z(i-1))); % calculate absorption

end

rad(i,j)=rad(i-1,j)*abs; %

end

end

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

% decide if and where to plot

ploc=find(plotix==n);  % is this one of the results we want to plot?

if not(isempty(ploc))  % then plot. “Ploc” is the location within all the plots

subplot(subr,subc,ploc),plot(v,rad(end,:))  % plot wavenumber against TOA emissive power

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

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

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

grid on

end

end

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

xlabel(‘pCO2 ppm’) % this x label is up for changing

====== end of RTE v0.2  =========

====== define_atmos ============

function [pr Tr rhor ztropo] = define_atmos(zr, Ts, ps)

% Calculate pressure, pr; temperature, Tr; density, rhor; from height, zr;

%  surface temperature, Ts; and surface pressure, ps

% By assuming a standard temperature profile to a fixed tropopause temp

% zr is a linear vector and starts at the surface

pr=zeros(1,length(zr)); % allocate space

Tr=zeros(1,length(zr)); % allocate space

rhor=zeros(1,length(zr)); % allocate space

% first calculate temperature profile

lapset=6.5/1000; % lapse rate in K/m

Ttropo=215; % temperature of tropopause

htropo=10000; % height of tropopause

Tstrato=270; % temperature of stratopause

zstrato=50000; % height of stratopause

ztropo=((Ts-Ttropo)/lapset); % height of bottom of tropopause in m

for i=1:length(zr)

if zr(i)<ztropo

Tr(i)=Ts-(zr(i)*lapset);  %  the troposphere

elseif (zr(i)>=ztropo && zr(i)<=(ztropo+htropo))

Tr(i)=Ttropo; % the tropopause

elseif (zr(i)>(ztropo+htropo) && zr(i)<zstrato)

Tr(i)=Ttropo+(zr(i)-(ztropo+htropo))*(Tstrato-Ttropo)/(zstrato-(ztropo+htropo)); % stratosphere

else

Tr(i)=Tstrato; % haven’t yet defined temperature above stratopause

% but this prevents an error condition

end

end

% pressure, p = ps * exp(-mg/R*integral(1/T)dz) from 0-z

mr=28.57e-3; % molar mass of air

R=8.31; % gas constant

g=9.8; % gravity constant

intzt=0; % sum the integral of dz/T in each iteration

dzr=zr(2)-zr(1); % dzr is linear so delta z is constant

pr(1)=ps; % surface condition

for i=2:length(zr)

intzt=intzt+(dzr/Tr(i));

pr(i)=ps*exp(-mr*g*intzt/R);

end

% density, rhor = mr.p/RT

rhor=mr.*pr./(R.*Tr);

end

======= end of define_atmos  =================

======= planckmv ==================

function ri = planckmv( v,t )

%  planck function for matrix of wavenumber, v (1/cm) and temp, t (K)

%  result in spectral emissive power (pi * intensity)

%  in W/(m^2.cm^-1)

%  from

http://pds-atmospheres.nmsu.edu/education_and_outreach/encyclopedia/planck_function.htm

[vv tt]= meshgrid(v,t);

ri =  3.7418e-8.*v.^3./(exp(vv.*1.4388./tt)-1);

end

============ end of planckmv =================

Read Full Post »

The Cool Skin of the Ocean

In the ensuing discussion on Does Back Radiation “Heat” the Ocean? – Part Four, the subject of the cool skin of the ocean surface came up a number of times.

It’s not a simple subject, but it’s an interesting one so I’m going to plough on with it anyway.

Introduction

The ocean surface is typically something like 0.1°C – 0.6°C cooler than the temperature just below the surface. And this “skin”, or ultra-thin region, is less than a 1mm thick.

Here’s a diagram I posted in the comments of Does Back Radiation “Heat” the Ocean? – Part Three:

Kawai & Wada (2007)

from Kawai & Wada (2007)

Figure 1

There is a lot of interest in this subject because of the question: “When we say ‘sea surface temperature’ what do we actually mean?“.

As many climate scientists note in their papers, the relevant sea surface temperature for heat transfer between ocean and atmosphere is the very surface, the skin temperature.

In figure 1 you can see that during the day the temperature increases up to the surface and then, in the skin layer, reduces again. Note that the vertical axis is a logarithmic scale.

Then at night the temperature below the skin layer is mostly all at the same temperature (isothermal). This is because the surface cools rapidly at night, and therefore becomes cooler than the water below, so sinks. This diurnal mixing can also be seen in some graphs I posted in the comments of Does Back Radiation “Heat” the Ocean? – Part Four.

Before we look at the causes, here are a series of detailed measurements from Near-surface ocean temperature by Ward (2006):

Figure 2

Note: The red text and arrow is mine, to draw attention to the lower skin temperature. The measurements on the right were taken just before midday “local solar time”. I.e., just before the sun was highest in the sky.

And in the measurements below I’ve made it a bit easier to pick out the skin temperature difference with blue text “Skin temp“. The blue value in each graph is what is identified as ΔTc in the schematic above. The time is shown as local solar time.

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

The measurements of the skin surface temperature were made by MAERI, a passive infrared radiometric interferometer. The accuracy of the derived SSTs from M-AERI is better than 0.05 K.

Below the skin, the high-resolution temperature measurements were measured by SkinDeEP, an autonomous vertical profiler. This includes the “sub-skin” measurement, from which the sea surface temperature was subtracted to calculate ΔTc (see figure 1).

The Theory

The existence of the temperature gradient is explained by the way heat is transferred: within the bulk waters, heat transfer occurs due to turbulence, but as the surface is approached, viscous forces dominate and molecular processes prevail. Because heat transfer by molecular conduction is less efficient than by turbulence, a strong temperature gradient is established across the boundary layer.

Ward & Minnett (2001)

Away from the interface the temperature gradient is quickly destroyed by turbulent mixing. Thus the cool-skin temperature change is confined to a region of thickness, which is referred to as the molecular sublayer.

Fairall et al (1996)

What do they mean?

Here’s an insight into what happens at fluid boundaries from an online textbook (thanks to Dan Hughes for letting me know about it) – this textbook is freely available online:

From "A Heat Transfer Textbook", by Prof Lienhard & Prof Lienhard (2008)

From "A Heat Transfer Textbook", by Prof Lienhard & Prof Lienhard (2008)

Figure 8

The idea behind turbulent mixing in fluids is that larger eddies “spawn” smaller eddies, which in turn spawn yet smaller eddies until you are up against an interface for that fluid (or until energy is dissipated by other effects).

In the atmosphere, for example, large scale turbulence moves energy across many 100’s of kilometers. A few tens of meters above the ground you might measure eddies of a few hundreds of meters in size, and in the last meter above the ground, eddies might be measured in cms or meters, if they exist at all. And by the time we measure the fluid flow 1mm from the ground there is almost no turbulence.

For some basic background over related terms, check out Heat Transfer Basics – Convection – Part One, with some examples of fluid flowing over flat plates, boundary layers, laminar flow and turbulent flow.

Therefore, very close to a boundary the turbulent effects effectively disappear, and heat transfer is carried out via conduction. Generally conduction is less effective than turbulence movement of fluids at heat transfer.

A Note on Very Basic Theory

The less effectively heat can move through a body, the higher the temperature differential needed to “drive” that heat through.

This is described by the equation for conductive heat transfer, which in (relatively) plain English says:

The heat flow in W/m² is proportional to the temperature difference across the body and the “conductivity” of the body, and is inversely proportional to the distance across the body

Now during the day a significant amount of heat moves up through the ocean to the surface. This is the solar radiation absorbed below the surface. Near the surface where turbulent mixing reduces in effectiveness we should expect to see a larger temperature gradient.

Taking the example of 1m down, if for some reason heat was not able to move effectively from 1m to the surface, then the absorbed solar radiation would keep heating the 1m depth and its temperature would keep rising. Eventually this temperature gradient would cause greater heat flow.

An example of a flawed model where heat was not able to move effectively was given in Does Back-Radiation “Heat” the Ocean? – Part Two:

A Flawed Model

Note how the 1m & 3m depth keep increasing in temperature. See that article for more explanation.

The Skin Layer in Detail

If the temperature increases closer to the surface, why does it “change direction” in the last millimeter?

In brief, the temperature generally rises in the last few meters as you get closer to the surface because hotter fluids rise. They rise because they are less dense.

So why doesn’t that continue to the very last micron?

The surface is where (almost) all of absorbed ocean energy is transferred to the atmosphere.

  • Radiation from the surface takes place from the top few microns.
  • Latent heat – evaporation of water into water vapor – is taken from the very top layer of the ocean.
  • Sensible heat is moved by conduction from the very surface into the atmosphere

And in general the ocean is moving heat into the atmosphere, rather than the reverse. The atmosphere is usually a few degrees cooler than the ocean surface.

Because turbulent motion is reduced the closer we get to the boundary with the atmosphere, this means that conduction is needed to transfer heat. This needs a temperature differential.

I could write it another way – because “needing a temperature differential” isn’t the same as “getting a temperature differential”.

If the heat flow up from below cannot get through to the surface, the energy will keep “piling up” and, therefore, keep increasing the temperature. Eventually the temperature will be high enough to “drive the heat” out to the surface.

The Simple 1-d Model

We saw a simple 1-d model in Does Back Radiation “Heat” the Ocean? – Part Four.

Just for the purposes of checking the theory relating to skin layers here is what I did to improve on it:

1. Increased the granularity of the model – with depths for each layer of: 100μm, 300μm, 1mm, 5mm, 20mm, 50mm, 200mm, 1m, 10m, 100m (note values are the lower edge of each layer).

2. Reduced the “turbulent conductivity” values as the surface was reached – instead of one “turbulent conductivity” value (used when the layer below was warmer than the layer above), these values were reduced closer to the surface, e.g. for the 100μm layer, kt=10; for the 300μm layer, kt=10; for the 1mm layer, kt=100; for the 5mm layer, kt=1000; for the 20mm layer, kt=100,000. Then the rest were 200,000 = 2×105 – the standard value used in the earlier models.

3. Reduced the time step to 5ms. This is necessary to make the model work and of course does reduce the length of run significantly.

The results for a 30 day run showed the beginnings of a cooler skin. And the starting temperatures for the top layer down to the 20mm layer were the same. The values of kt were not “tuned” to make the model work, I just threw some values in to see what happened.

As a side note for those following the discussion from Part Four, the ocean temperature also increased for DLR increases with these changes.

Now I can run it for longer but the real issue is that the model is not anywhere near complex enough.

Further Reading on Complexity

There are some papers for people who want to follow this subject further. This is not a “literature review”, just some papers I found on the journey. The subject is not simple.

Free

Saunders, Peter M. (1967), The Temperature at the Ocean-Air InterfaceJ. Atmos. Sci.

Tu and Tsuang (2005), Cool-skin simulation by a one-column ocean model, Geophys. Res. Letters

Paywall

McAlister, E. D., and W. McLeish (1969), Heat Transfer in the Top Millimeter of the OceanJ. Geophys. Res.

Fairall et al, reference below

GA Wick, WJ Emery, LH Kantha & P Schlussel (1996), The behavior of the bulk-skin sea surface temperature difference under varying wind speed and heat fluxJournal of Physical Oceanography

Hartmut Grassl, (1976), The dependence of the measured cool skin of the ocean on wind stress and total heat flux, Boundary Layer Meteorology

Conclusion

The temperature profile of the top mm of the ocean is a challenging subject. Tu & Tsuang say:

Generally speaking, the structure of the viscous layer is known to be related to the molecular viscosity, surface winds, and air-sea flux exchanges. Both Saunders’ formulation [Saunders, 1967; Grassl, 1976; Fairall et al.,1996] and the renewal theory [Liu et al., 1979; Wick et al.,1996; Castro et al., 2003; Horrocks et al., 2003] have been developed and applied to study the cool-skin effect.

But the exact factors and processes determining the structure is still not well known.

However, despite the complexity, an understanding of the basics helps to give some insight into why the temperature profile is like it is.

I welcome commenters who can make the subject easier to understand. And also commenters who can explain the more complex elements of this subject.

References

A Heat Transfer Textbook, by Prof Lienhard & Prof Lienhard, Phlogiston Press, 3rd edition (2008)

Cool-skin and warm-layer effects on sea surface temperature, Fairall, Bradley, Godfrey, Wick, Edson & Young, Journal of Geophysical Research (1996)

Near-surface ocean temperature, Ward, Journal of Geophysical Research (2006)

An Autonomous Profiler for Near Surface Temperature Measurements, Ward & Minnett, Accepted for the Proceedings Gas Transfer at Water Surfaces 4th International Symposium (2000)

Read Full Post »

In Part One we saw how the ocean absorbed different wavelengths of radiation:

  • 50% of solar radiation is absorbed in the first meter, and 80% within 10 meters
  • 50% of “back radiation” (atmospheric radiation) is absorbed  in the first few microns (μm).

This is because absorption is a strong function of wavelength and atmospheric radiation is centered around 10μm, while solar radiation is centered around 0.5μm.

In Part Two we considered what would happen if back radiation only caused evaporation and removal of energy from the ocean surface via the latent heat. The ocean surface would become much colder than it obviously is. That is a very simple “first law of thermodynamics” problem. Then we looked at another model with only conductive heat transfer between different “layers” in the ocean. This caused various levels below the surface to heat to unphysical values. It is clear that turbulent heat transport takes place from lower in the ocean. Solar energy reaches down many meters heating the ocean from within – hotter water expands and so rises – moving heat by convection.

In Part Three we reviewed various experimental results showing how the temperature profile (vs depth) changes during the diurnal cycle (day-night-day) and with wind speed. This demonstrates very clearly how much mixing goes on in the ocean.

The Different Theories

This series of articles was inspired by the many people who think that increases in back radiation from the atmosphere will have no effect (or an unnoticeable effect) on the temperature of the ocean depths.

So far, no evidence has so far been brought forward for the idea that back radiation can’t “heat” the ocean (see note 1 at the end), other than the “it’s obvious” evidence. At least, I am unaware of any stronger arguments. Hopefully as a result of this article advocates can put forward their ideas in more detail in response.

I’ll summarize the different theories as I’ve understood them. Apologies to anyone who feels misrepresented – it’s quite possible I just haven’t heard your particular theory or your excellent way of explaining it.

Hypothesis A – Because the atmospheric radiation is completely absorbed in the first few microns it will cause evaporation of the surface layer, which takes away the energy from the back radiation as latent heat into the atmosphere. Therefore, more back-radiation will have zero effect on the ocean temperature.

Hypothesis B – Because the atmospheric radiation is completely absorbed in the first few microns it will be immediately radiated or convected back out to the atmosphere. Heat can’t flow downwards due to the buoyancy of hotter water. Therefore, if an increase in back radiation occurs (perhaps due to increases in inappropriately-named “greenhouse” gases) it will not “heat” the ocean = increase the temperature of the ocean below the surface.

For other, more basic objections about back radiation, see Note 2 (at the end).

I believe that Part Two showed that Hypothesis A was flawed.

I would like to propose a different hypothesis:

Hypothesis C – Heat transfer is driven by temperature differences. For example, conduction of heat is proportional to the temperature difference across the body that the heat is conducted through.

Solar radiation is absorbed from the surface through many meters of the ocean. This heats the ocean below the surface which causes “natural convection” – heated bodies expand and therefore rise. So solar energy has a tendency to be moved back to the surface (this was demonstrated in Part Two).

The more the surface temperature increases, the less temperature difference there will be to drive this natural convection. And, therefore, increases in surface temperature can affect the amount of heat stored in the ocean.

Clarification from St.Google: HypothesisA supposition or proposed explanation made on the basis of limited evidence as a starting point for further investigation

An Excellent Question

In Part Three, one commenter asked an excellent question:

Some questions from an interested amateur.
Back radiation causes more immediate evaporation and quicker reemission of LWR than does a similar amount of solar radiation.

Does that mean that the earth’s temperature should be more sensitive to a given solar forcing than it would be to an equal CO2 forcing?

What percentage CO2 forcing transfers energy to the oceans compared to space and the atmosphere?

How does this compare with solar forcing?

Is there a difference between the effect of the sun and the back radiation when they are of equal magnitude? This, of course, pre-supposes that Hypothesis C is correct and that back radiation has any effect at all on the temperature of the ocean below the surface.

So the point is this - even if Hypothesis C is correct, there may still be a difference between the response of the ocean temperatures below the surface – for back radiation compared with solar radiation.

So I set out to try and evaluate these two questions:

  1. Can increases in back radiation affect the temperature of the ocean below the surface? I.e., is Hypothesis C supported against B?
  2. For a given amount of energy, is there a difference between solar forcing and back radiation forcing?

And my approach was to use a model:

Oh no, a model! Clearly wrong then, and a result that can’t fool anyone..

For a bit of background generally on models, take a look at the Introduction in Models On – and Off – the Catwalk.

Here is one way to think about a model

The idea of a model is to carry out some calculations when doing them in your head is too difficult

A model helps us see the world a bit more clearly. At these point I’m not claiming anything other than they help us see the effect of the well-known and undisputed laws of heat transfer on the ocean a little bit more clearly.

Ocean Model

The ocean model under consideration is about a billion times less complex than a GCM. It is a 1-d model with heat flows by radiation, conduction and, in a very limited form, convection.

Here is a schematic of the model. I thought it would be good to show the layers to scale but that means the thicker layers can’t be shown (not without taking up a ridiculous amount of blank screen space) – so the full model, to scale, is 100x deeper than this:

Figure 1

To clarify – the top layer is at temperature, T1, the second layer at T2, even though these values aren’t shown.

The red arrows show conducted or convected heat. They could be in either direction, but the upwards is positive (just as a convention). Obviously, only a few of these are shown in the schematic – there is a heat flux between each layer.

1. Solar and back radiation are modeled as sine waves with the peak at midday. See the graph “Solar and Back Radiation” in Part Two for an example.

2. Convected heat is modeled with a simple formula:

H=h(T1-Tair), where Tair = air temperature, T1 = “surface” temperature, h = convection coefficient = 25 W/m².K.

Convected heat can be in either direction, depending on the surface and air temperature. The air temperature is assumed constant at 300K, something we will return to.

3. Radiation from the surface:

E = εσT4 – the well-known Stefan-Boltzmann equation, and ε = emissivity

For the purposes of this simple model ε = 1. So is absorptivity for back radiation, and for solar radiation. More on these assumptions later.

4. Heat flux between layers (e.g. H54 in the schematic) is calculated using the temperature values for the previous time step for the two adjacent layers then using the conducted heat formula: q” = k.(T5-T4)/d54, where k= conductivity, and d54 = distance between center of each layer 5 to the center of layer 4.

For still water, k = 0.6 W/m.K – a very low value as water is a poor conductor of heat.

In this model at the end of each time step, the program checks the temperature of each layer. If T5 > T4 for example, then the conductivity between these layers for the next time step is set to a much higher value to simulate convection. I used a value for stirred water that I found in a textbook: kt = 2 x 105 W/m.K. What actually happens in practice is the hotter water rises taking the heat with it (convection). Using a high value of conductivity produces a similar result without any actual water motion.

For interest I did try lower values like 2 x 10³ W/m.K and the 1m layer, for example, ended up at a higher temperature than the layers above it. See the more detailed explanation in Part Two.

5. In Part Three I showed results from a number of field experiments which demonstrated that the ocean experiences mixing due to surface cooling at night, and due to high winds. The mixing due to surface cooling is automatically taken account of in this model (and we can see it in the results), but the mixing due to the winds “stirring” the ocean is not included. So we can consider the model as being “under light winds”. If we had a model which evaluated stronger winds it would only make any specific effects of back radiation less noticeable. So this is the “worst case” – or the “highlighting back radiation’s special nature” model.

Problems of Modeling

Some people will already know about the many issues with numerical models. A very common one is resolving small distances and short timescales.

If we want to know the result over many years we don’t really want to have the iterate the model through time steps of fractions of a second. In this model I do have to use very small time steps because the distance scales being considered range from extremely small to quite large – the ocean is divided into thin slabs of 5mm, 15mm.. through to a 70m slab.

If I use a time step which is too long then too much heat gets transferred from the layers below the surface to the 5mm surface layer in the one time step, the model starts oscillating - and finally “loses the plot”. This is easy to see, but painful to deal with.

But I thought it might be interesting for people to see the results of the model over five days with different time steps. Instead of having the model totally “lose the plot” (=surface temperature goes to infinity), I put a cap on the amount of heat that could move in each time step for the purposes of this demonstration.

You can see four results with these time steps (tstep = time step, is marked on the top left of each graph):

  • 3 secs
  • 1 sec
  • 0.2 sec
  • 0.05 sec

Figure 2 – Click for a larger image

I played around with many other variables in the model to see what problems they caused..

The Tools

The model is written in Matlab and runs on a normal PC (Dell Vostro 1320 laptop).

To begin with there were 5 layers in the model (values are depth from the surface to the bottom edge of each layer):

  • 5 mm
  • 50 mm
  • 1 m
  • 10 m
  • 100 m

I ran this with a time step of 0.2 secs and ended up doing up to 15-year runs.

In the model runs I wanted to ensure that I had found a steady-state value, and also that the model conserved energy (first law of thermodynamics) once steady state was reached. So the model included a number of “house-keeping” tests so I could satisfy myself that the model didn’t have any obvious errors and that equilibrium temperatures were reached for each layer.

For 15 year runs, 5 layers and 0.2s time step the run would take about two and a half hours on the laptop.

I find that quite amazing – showing how good Matlab is. There are 31 million seconds in a year, so 15 years at 0.2 secs per step = 2.4 billion iterations. And each iteration involves looking up the solar and DLR value, calculating 7 heat flow calculations and 5 new temperatures. All in a couple of hours on a laptop.

Well, as we will see, because of the results I got I thought I would check for any changes if there were more layers in my model. So that’s why the 9-layer model (see the first diagram) was created. For this model I need an even shorter time step – 0.1 secs and so long model runs start to get painfully long..

Results

Case 1: The standard case was a peak solar radiation, S, of 600 W/m² and back radiation, DLR of 340 with a 50 W/m² variation day to night (i.e., max of 390 W/m², min of 290 W/m²).

Case 2a: Add 10 W/m² to the peak solar radiation, keep DLR the same. Case 2b – Add 31.41 W/m² to solar.

Case 3a: Keep solar radiation the same, add 3.14 W/m² to DLR. This is an equivalent amount of energy per day to case 2, see note 3. Case 3b – Add 10 W/m² to DLR.

Many people are probably asking, “Why isn’t case 3a – Add 10 W/m² to DLR?”

Solar radiation only occurs for 12 out of the 24 hours, while DLR occurs 24 hours of the day. And the solar value is the peak, while the DLR value is the average. It is a mathematical reason explained further in Note 3.

The important point is that for total energy absorbed in a day, case 2a and 3a are the same, and case 2b and 3b are the same.

Let’s compare the average daily temperature in the top layer, 1m, 10m and 100m layer for the three cases (note: depths are from the surface to the bottom of each layer; and only 4 layers of the 5 were recorded):

Figure 3

The time step (tstep) = 0.2s.

The starting temperatures for each layer were the same in all cases.

Now because the 4 year runs recorded almost identical values for solar vs DLR forcing, and because the results had not quite stabilized, I then did the 15 year run and also recorded the temperature to the 4 decimal places shown. This isn’t because the results are this accurate – this is to see what differences, if any, exist between the two different scenarios.

The important results are:

  1. DLR increases cause temperature increases at all levels in the ocean
  2. Equivalent amounts of daily energy into the ocean from solar and DLR cause almost exactly the same temperature increase at each level of the ocean – even though the DLR is absorbed in the first few microns and the solar energy in the first few meters
  3. The slight difference in temperature may be a result of “real physics” or may be an artifact of the model

And perhaps 5 layers is not enough?

Therefore, I generated the 9-layer model, as shown in the first diagram in this article. The 15-year model runs on the 9-layer model produced these results:

Figure 4

The general results are similar to the 5-layer model.

The temperature changes have clearly stabilized, as the heat unaccounted for (inputs – outputs) on the last day = 41 J/m². Note that this is Joules, not Watts, and is over a 24 hour period. This small “unaccounted” heat is going into temperature increases of the top 100m of the ocean. (“Inputs – outputs” includes the heat being transferred from the model layers down into the ocean depths below 100m).

If we examine the difference in temperature for the bottom 30-100m deep level for case 2b vs 3b we see that the temperature difference after 15 years = 0.011°C. For a 70m thick layer, this equates to an energy difference = 3.2 x 106 J, which, over 15 years, = 591 J/m².day = 0.0068 W/m². This is spectacularly tiny. It might be a model issue, or it might be a real “physics difference”.

In any case, the model has demonstrated that DLR increases vs solar increases cause almost exactly the same temperature changes in each layer being considered.

For interest here are the last 5 days of the model (average hourly temperatures for each level) for case 3b:

Figure 5

and for case 2b:

Figure 6

Pretty similar..

Results – Convection and Air Temperature

In the model results up until now the air temperature has been at 300K (27°C) and the surface temperature of the ocean has been only a few degrees higher.

The model doesn’t attempt to change the air temperature. And in the real world the atmosphere at the ocean surface and the surface temperature are usually within a few degrees.

But what happens in our model if real world situations cool the ocean surface more? For example, higher temperatures locally create large convective currents of rising hot air which “sucks in” cooler air from another area.

What would be the result? A higher “instantaneous” surface temperature from higher back radiation might be “swept away” into the atmosphere and “lost” from the model.. This might create a different final answer for back radiation compared with solar radiation.

It seemed to be worth checking out, so I reduced the air temperature to 285K (from 300K) and ran the model for one year from the original starting temperatures (just over 300K). The result was that the ocean temperature dropped significantly, demonstrating how closely the ocean surface and the atmosphere (at the ocean surface) are coupled.

Using the end of the first year as a starting temperature, I ran the model for 5 years for case 1, 2a and 3a (each with the same starting temperature):

Figure 7

Once again we see that back radiation increases do change the temperatures of the ocean depths – and at almost identical values to the solar radiation changes.

Here is a set of graphs for one of the 5-year model runs for this lower air temperature, also demonstrating how the lower air temperature pulls down the ocean surface temperature:

Figure 8 – Click for a larger image

The first graph shows how the average daily temperature changes over the full time period – making it easy to see equilibrium being reached. The second graph shows the hourly average temperature change for the last 5 days. The last graph shows the heat which is either absorbed or released within the ocean in temperature changes. As zero is reached it means the ocean is not heating up or cooling down.

Inaccuracies in the Model

We can write a lot on the all the inaccuracies in the model. It’s a very rudimentary model. In the real world the hotter tropical / sub-tropical oceans transfer heat to higher latitudes and to the poles. So does the atmosphere. A 1-d model is very unrealistic.

The emissivity and absorptivity of the ocean are set to 1, there are no ocean currents, the atmosphere doesn’t heat up and cool down with the ocean surface, the solar radiation value doesn’t change through the year, the top layer was 5mm not 1μm, the cooler skin layer was not modeled, a number of isothermal layers is unphysical compared with the real ocean of continuously varying temperatures..

However, what a nice simple model tells us is how energy only absorbed in the top few microns of the ocean can affect the temperature of the ocean much lower down.

It’s obvious“, I could say.

Conclusion

My model could be wrong – for example, just a mistake which means it doesn’t operate how I have described it. The many simplifications of the model might hide some real world physics effect which means that Hypothesis C is actually less likely than Hypothesis B.

However, if the model doesn’t contain mistakes, at least I have provided more support for Hypothesis C – that the back radiation absorbed in the very surface of the ocean can change the temperature of the ocean below, and demonstrated that Hypothesis B is less likely.

I look forward to advocates of Hypothesis B putting forward their best arguments.

Update - Code files saved here

Notes

Note 1 – To avoid upsetting the purists, when we say “does back-radiation heat the ocean?” what we mean is, “does back-radiation affect the temperature of the ocean?”

Some people get upset if we use the term heat, and object that heat is the net of the two way process of energy exchange. It’s not too important for most of us. I only mention it to make it clear that if the colder atmosphere transfers energy to the ocean then more energy goes in the reverse direction.

It is a dull point.

Note 2 - Some people think that back radiation can’t occur at all, and others think that it can’t affect the temperature of the surface for reasons that are a confused mangle of the second law of thermodynamics. See Science Roads Less Travelled and especially Amazing Things we Find in Textbooks – The Real Second Law The Real Second Law of Thermodynamics and The Three Body Problem. And for real measurements of back radiation, see The Amazing Case of “Back Radiation” -Part One.

Note 3 – If we change the peak solar radiation from 600 to 610, this is the peak value and only provides an increase for 12 out of 24 hours. By contrast, back radiation is a 24 hour a day value. How much do we have to change the average DLR value to provide an equivalent amount of energy over 24 hours?

If we integrate the solar radiation for the before and after cases we find the relationship between the value for the peak of the solar radiation and the average of the back radiation = π (3.14159). So if the DLR increase = 10, the peak solar increase to match = 10 x π = 31.4159; and if the solar peak increase = 10, the DLR increase to match = 10/π = 3.1831.

If anyone would like this demonstrated further please ask and I will update in the comments. I’m sure I could have made this easier to understand than I actually have (haven’t).

Read Full Post »

A while ago we looked at some basics in Heat Transfer Basics – Part Zero.

Equations aren’t popular but a few were included.

As a recap, there are three main mechanisms of heat transfer:

  • conduction
  • convection
  • radiation

In the climate system, conduction is generally negligible because gases and liquids like water don’t conduct heat well at all. (See note 2).

Convection is the transfer of heat by bulk motion of a fluid. Motion of fluids is very complex, which makes convection a difficult subject.

If the motion of the fluid arises from an external agent, for instance, a fan, a blower, the wind, or the motion of a heated object itself, which imparts the pressure to drive the flow, the process is termed forced convection.

If, on the other hand, no such externally induced flow exists and the flow arises “naturally” from the effect of a density difference, resulting from a temperature or concentration difference in a body force field such as gravity, the process is termed natural convection. The density difference gives rise to buoyancy forces due to which the flow is generated..

The main difference between natural and forced convection lies in the mechanism by which flow is generated.

From Heat Transfer Handbook: Volume 1, by Bejan & Kraus (2003).

The Boundary Layer

The first key to understanding heat transfer by convection is the boundary layer. A typical example is a fluid (e.g. air, water) forced over a flat plate:

From Incropera & DeWitt (2007)

From Incropera & DeWitt (2007)

This first graphic shows the velocity of the fluid. The parameter u is the velocity (u) at infinity (∞) – or in layman’s terms, the velocity of the fluid “a long way” from the surface of the plate.

Another way to think about u – it is the free flowing fluid velocity before the fluid comes into contact with the plate.

Take a look at the velocity profile:

At the plate the velocity is zero. This is because the fluid particles make contact with the surface. In the “next layer” the particles are slowed up by the boundary layer particles. As you go further and further out this effect of the stationary plate is more and more reduced, until finally there is no slowing down of the fluid.

The thick black curve, δ, is the boundary layer thickness. In practice this is usually taken to be the point where the velocity is 99% of its free flowing value. You can see that just at the point where the fluid starts to flow over the plate – the boundary layer is zero. Then the plate starts to slow the fluid down and so progressively the boundary layer thickens.

Here is the resulting temperature profile:

From Incropera & DeWitt (2007)

From Incropera & DeWitt (2007)

In this graphic T is the temperature of the “free flowing fluid” and Ts is the temperature of the flat plat which (in this case) is higher than the free flowing fluid temperature. Therefore, heat will transfer from the plate to the fluid.

The thermal boundary layer, δt, is defined in a similar way to the velocity boundary layer, but using temperature instead.

How does heat transfer from the plate to the fluid? At the surface the velocity of the fluid is zero and so there is no fluid motion.

At the surface, energy transfer only takes place by conduction (note 1).

In some cases we also expect to see mass transfer – for example, air over a water surface where water evaporates and water vapor gets carried away. (But not with air over a steel plate).

From Incropera & DeWitt (2007)

From Incropera & DeWitt (2007)

So a concentration boundary layer develops.

Newton’s Law of Cooling

Many people have come across this equation:

q” = h(Ts – T)

where q” = heat flux in W/m², h is the convection coefficient, and the two temperatures were defined above

The problem is determining the value of h.

It depends on a number of fluid properties:

  • density
  • viscosity
  • thermal conductivity
  • specific heat capacity

But also on:

  • surface geometry
  • flow conditions

Turbulence

The earlier examples showed laminar flow. However, turbulent flow often develops:

Flow in the turbulent region is chaotic and characterized by random, three-dimensional motion of relatively large parcels of fluid.

Check out this very short video showing the transition from laminar to turbulent flow.

What determines whether flow is laminar or turbulent and how does flow become turbulent?

The transition from laminar to turbulent flow is ultimately due to triggering mechanisms, such as the interaction of unsteady flow structures that develop naturally within the fluid or small disturbances that exist within many typical boundary layers. These disturbances may originate from fluctuations in the free stream, or they may be induced by surface roughness or minute surface vibrations

from Incropera & DeWitt (2007).

Imagine treacle (=molasses) flowing over a plate. It’s hard to picture the flow becoming turbulent. That’s because treacle is very viscous. Viscosity is a measure of how much resistance there is to different speeds within the fluid – how much “internal resistance”.

Now picture water moving very slowly over a plate. Again it’s hard to picture the flow becoming turbulent. The reason in this case is because inertial forces are low. Inertial force is the force applied on other parts of the fluid by virtue of the fluid motion.

The higher the inertial forces the more likely fluid flow is to become turbulent. The higher the viscosity of the fluid the less likely the fluid flow is to become turbulent – because this viscosity damps out the random motion.

The ratio between the two is the important parameter. This is known as the Reynolds number.

Re = ρux / μ

where ρ = density, u = free stream velocity, x is the distance from the leading edge of the surface and μ = dynamic viscosity

Once Re goes above around 5 x 105 (500,000) flow becomes turbulent.

For air at 15°C and sea level, ρ=1.2kg/m³ and μ=1.8 x 10-5 kg/m.s

Solving this equation for these conditions, gives a threshold value of ux > 7.5 for turbulence.. This means that if the wind speed (in m/s) x the length of surface over which the wind flows (in m) is greater than 7.5 we will get turbulent flow.

For example, a slow wind speed of 1 m/s (2.2 miles / hour) over 7.5 meters of surface will produce turbulent flow. When you consider the wind blowing over many miles of open ocean you can see that the air flow will almost always be turbulent.

The great physicist and Nobel Laureate Richard Feynman called turbulence the most important unsolved problem of classical physics.

In a nutshell, it’s a little tricky. So how do we determine convection coefficients?

Empirical Measurements & Dimensionless Ratios

Calculation of the convection heat transfer coefficient, h, in the equation we saw earlier can only be done empirically. This means measurement.

However, there are a whole suite of similarity parameters which allow results from one situation to be used in “similar circumstances”.

It’s not an easy subject to understand “intuitively” because the demonstration of these similarity parameters (e.g., Reynolds, Prandtl, Nusselt and Sherwood numbers) relies upon first seeing the differential equations governing fluid flow and heat & mass transfer – and then the transformation of these equations into a dimensionless form.

As the simplest example, the Reynolds number tells us when flow becomes turbulent regardless of whether we are considering air, water or treacle.

And a result for one geometry can be re-used in a different scenario with similar geometries.

Therefore, many tables and standard empirical equations exist for standard geometries – e.g. fluid flow over cylinders, banks of pipes.

Here are some results for air flow over a flat isothermal plate (isothermal = all at the same temperature) – calculated using empirically-derived equations:

Click for a larger view

The 1st graph shows that the critical Reynolds number of 5×105 is reached at 1.3m. The 2nd graph shows how the boundary layer grows under first laminar flow, then second under turbulent flow – see how it jumps up as turbulent flow starts. The 4th graph shows the local convection coefficient as a function of distance from the leading edge – as well as the average value across the 2m of flat plate.

Conclusion

Not much of a conclusion yet, but this article is already long enough. In the next article we will look at the experimental results of heat transfer from the ocean to the atmosphere.

Notes

Note 1 – Heat transfer by radiation might also take place depending on the materials in question.

Note 2 – Of course, as explained in the detailed section on convection, heat cannot be transferred across a boundary between a surface and a fluid by convection. Conduction is therefore important at the boundary between the earth’s surface and atmosphere.

Read Full Post »

Follow

Get every new post delivered to your Inbox.

Join 289 other followers