Saturday 27 December 2014

Calculating insolation as a function of date and latitude.

In preparation for some work I will be doing I've been working on calculation of insolation as a function of time and latitude. In this post I will explain the method and compare results to a readily obtainable plot of insolation at various latitudes.

Note, all of the following calculations and the comparison with another source are for top of atmosphere and do not take into account absorption/scattering by the atmosphere, or reflection from clouds.




Insolation is a term for the solar radiation received by the earth, most of the insolation is in the visible spectrum, indeed that is the reason our eyes see in that spectrum (hence visible spectrum). The calculations have been carried out in an Excel 2007 spreadsheet, this should work in earlier and later versions of Excel, whether it will work in other programs (such as Libre Office) I do not know. The spreadsheet is available here. The maths and theory behind the calculation has been obtained from these lecture notes. Thanks are due once again to an academic making their course material available to the public!

Before going into detail, here is a screenshot of the upper left corner of the spreadsheet. Note that the statement 'assume circular orbit' is something I forgot to remove from the initial version of the calculation, an elliptical orbit is assumed, as will be explained.


I've colour coded it so as to try to make the various elements more clear. The largest single region, ranging from greens to reds is the actual results of calculated insolation. To the left of that is a pink region with latitude in degrees and in radians.

The constants used are at the top of the sheet, tilt, orbital period, eccentricity of the earth, also solar constant (1361W/m^2) and vernal equinox (20 March).

The year day is given in five day increments, it's easy to change that if needed because of the way the formula has been implemented. Above that is the ratio of the average earth to sun distance (ro) and the actual distance factoring in the eccentricity of the orbit of the earth. The vernal day is the day of a year centred on the spring equinox (20 March ~ Day 80). And the declination gives the variation of the angle of the sun throughout the year (in radians) due to the tilt of the earth.

Now for the nuts and bolts.

The length of day is H/Pi (I don't have the easy facility to use greek notation here, the English version of such notation is italicised from here on). H itself is the hour angle at sunset. H is calculated using latitude (phi) and the declination of the sun (delta), presented in terms of Excel formulae the formula for H is:

H = ACOS(-TAN(phi)*TAN(delta))

And here we hit a problem, in the Arctic the sun never rises for periods in the winter, and never sets for periods in the summer, the result is that the -TAN(phi)*TAN(delta) produces unphysical results which are greater than +/-1 so the ACOS (inverse cosine) generates a numeric error at high latitudes in summer and winter. A cosine always returns a result between +1 and -1 so the inverse cosine cannot be fed with a number outside that range.

H = ACOS(
IF(ABS(-1*TAN(phi)*TAN(delta))>1,SIGN(-1*TAN(phi)*TAN(delta)),-1*TAN(phi)*TAN(delta))
)

Already that looks a bit daunting, but it isn't too bad really. ABS() is a function that produces the absolute numeric value, without the sign of the number. SIGN() is a function that produces a +1 if the sign of the number is positive, and a -1 if the sign of the number is negative. So the above is described by:

If ABS(-TAN(phi)*TAN(delta)) is greater than 1 then use the sign of the "-TAN(phi)*TAN(delta)" part, otherwise use the actual value "-TAN(phi)*TAN(delta)".

This means the result is always equal to or less than 1 and stops the numeric errors at the polar regions.

From the lecture notes linked to above the following equation is provided, it describes the effect of the elliptical orbit of the earth on the earth to sun distance, which obviously impacts received insolation.


This gives the basics of a term (ro/r)^2 in the main equation, described below. Using the latter of the above two equations I use the following for the (ro/r)^2 term:

(ro/r)^2 = 1/(1-Eccentricity*COS(2*PI()*((YearDay-4)/OrbitalPeriod)))

Note that the average distance between earth and sun (ro) and the distance adjusted for the ellipse of the earth's orbit (r) are not actually used, because both 'r' terms cancel out. The result is the ratio of the two which is all we need.

So from the lecture notes linked to at the top of this post we have the main equation.


We have the term 'So' (Solar Constant) as the insolation from the sun received in space at the disantance from earth to sun perpendicular to the sun. 1361W/m^2 is used for this, sourced from the Wikipedia page on Solar Constant. We have the (ro/r)^2 term, which is implemented as a row in the spreadsheet (row 7), term H is as described above and is implemented within the equation in each spreadsheet cell where insolation is calculated. Latitude (phi) is calculated in radians and provided in column E of the spreadsheet. Declination (delta) is calculated and provided in row 10 of the spreadsheet.

The variable setting rows and columns then provide the drivers for a rectangular matrix of calculation of insolation as a function of latitude (vertical) and time (horizontal). The equation has been set up so that copying and pasting, also fill down, work even if the size of the range is changed to increase or decrease time or latitudinal resolution (using the $ operator in cell notation).

The main equation as used in the sheet is:

=(So*ro/r^2/PI())*
(
ACOS(IF(ABS(-1*TAN(delta)*TAN(phi))>1,SIGN(-1*TAN(delta)*TAN(phi)),-1*TAN(delta)*TAN(phi)))
*SIN(phi)*SIN(delta)+
SIN(ACOS(IF(ABS(-1*TAN(delta)*TAN(phi))>1,SIGN(-1*TAN(delta)*TAN(phi)),-1*TAN(delta)*TAN(phi))))
*COS(phi)*COS(delta)
)

The 'H' term being highlighted in blue to aid reading.

So how does it all work out? A commonly used presentation of insolation as a function of time and latitude is available from PhysicalGeography.net. This plot uses a set of latitude bands.


Using that as a reference of the correct answer worked out by a proper scientist (I presume), I can compare it with my own amateur calculation.


The results are pretty close, but not exactly (is anything ever so?). Note that for 60 and 30degN in my calculation the insolation is the same in mid summer, whereas in the official plot, 60degN is slightly lower. Also the mid summer peak for 90degN is too high by a small margin in my plot. However overall the agreement is good, and for my purposes what looks like a few percent of error at high latitudes is not likely to be a problem.

28 comments:

Anonymous said...

These are incoming solar flux values and not total energy received, right? Flux in Watts per metre squared and energy in Watt hours. Do we multiply the above values by 24 hours or by day length for the given period?

Chris Reynolds said...

They're watts per square metre. Never used watt hours though, you'd need to look that up.

The total solar energy is multiplied by factors including day length (H). So I think the correct interpretation is average over a 24 hour period. For joules that's w/squre metre x 3600 x 24.

Actually watt hour would probably be watts x3600?

Chris Reynolds said...

Just checked, Watt-hour is actually Watts X Hours. Not Watts X 3600

Mack said...

I think if you average these graphs out it comes to 342w/sq.m average per year over the whole globe, doesn't it Chris? ie solar insolation over the Earth's SURFACE is 342w/sq.m. This does not sit well with Trenberth's Earth energy budget diagrams where we have 161 or 163 or 168w/sq.m (seen various numbers in that range)solar radiation striking the Earth's surface.
Methinks you are right Chris, and Trenberth and all the rest with their looney Earth Energy Budget cartoons have effed up big time.
What say you?

Chris Reynolds said...

Mack,

The average insolation for my total table is just over 295W/m^2, a reduced table of 10 day increments I've been using for another blog post is just under 295W/m^2. In other words, about 295W/m^2, not far off the Kiehl/Trenberth assumption of 340W/m^2 insolation:
http://scied.ucar.edu/sites/default/files/users/lisagard/radiation_budget_kiehl_trenberth_2008_big.jpg

The calculations in this blog post are for top of atmosphere, not allowing for clouds or scattering in the atmosphere. I've not included those because that looks far too difficult for me.

Anyway Kiehl/Trenberth estimate about 79W/m^2 reflected by clouds and atmosphere globally, and about 78W/m^2 absorbed by the atmosphere. Call those 80 each, totalling 160W/m^2

295 - 160 = 135W/m^2

That's in the same ball park as the Kiehl/Trenberth figure of 161W/m^2.

The source of the confusion here is myself. I omitted to explicitly say what the insolation being calculated is, i.e. it's not surface, it's top of atmosphere. I assumed it was obvious, when it clearly wasn't.

I've put a note making that clear at the start of the post.

Sorry for the confusion.

Mack said...

"Sorry for the confusion"
There's no need to apologise for that little bit of confusion Chris,..but what you do need to ask is why somebody like me should come along and create this little bit of confusion to a very bright person like yourself, who obviously has a good handle on maths and physics..but somewhere in all those lectures,tutorials and educational materials you've assimilated, the fact you have to say that these watts/sq.m are landing on the TOA and not the Earth's surface,should suddenly jolt you. I mean, look at your link to PhysicalGeography.net .Nowhere there is there any mention of the TOA. In fact it looks extremely like they talk about the Earth's surface and watts/sq m. impinging thereon. So what's going on?
Well what's going on is great confusion and coverup by every learning institute as to which surface they are actually referring to with respect to the "greenhouse" theory.
I'm a denialist Chris...so watch out ! :)

Chris Reynolds said...

Mack,

I really think 'cover up' is too strong. In a complex subject it is all too easy to forget that the public have an interest but don't have the background learning to see what might seem obvious to a few. My error here was exactly that. Having read the source document it was obvious to me that nowhere was atmospheric absorption/scattering addressed, therefore the results were for top of atmosphere.

But not everyone reads all the post. I have done the same with papers where the maths is over my head, I read the abstract, introduction, conclusions and look at the graphs. It can be a risky thing to do and never means I really understand the paper concerned.

Your attitude to (anthropogenic) GW science was clear. But mainly you made a reasonable point, so I addressed it and didn't delete. Note that I am equally as hard on the extreme alarmists who lose touch with reality. I can be rude and impatient sometimes, but I do dish it out to both sides. :)

Mack said...

Chris, Thankyou for addressing my reasonable point and not deleting it. I can tell you have a scientific, truth-seeking and open mind....still !? :) and on the strength of that I wish to prevale upon you to investigate this link...
http://www.principia-scientific.org/the-greenhouse-effect-and-the-infrared-radiative-structure-of-the-earth-s-atmosphere.html#comment-9486
where I made only 2 comments just the other day to a Mr Doug Cotton, whom you may have heard of, towards the end of this thread. If you were to follow up my links, you'll discover a very interesting read which should appeal to you. (looks like Doug and this guy geran are still fighting it out, so leave them to it).

Chris Reynolds said...

Mack,

Up until 2007 I was rather sceptical of climate change and the human cause, in 2007 I read a paper which pretty much blew my solar-GW link out of the water (Wild et al "Impact of global dimming and brightening on global warming"). I've skimmed the contents on that page and see nothing that impresses me.

There are three possibilities:

1) That the people there have got the physics right and the rest of the physicists are wrong.

2) That the people there have got the physics right and the rest of the physicists are engaging in a mass deception.

3) The people saying they've shown physics wrong just don't really get the physics that persuades the vast majority of physicists.

1 is unlikely because as a reader of climate science (now mainly the Arctic) I can assure you that there are very intelligent people working on the problem.

2 is inherently implausible because of the competitiveness of scientists. Also because most people are honest and scientists comprise a broad spectrum of human behavioural types.

3 is most likely. From my own experience, there have been but one or two occasions in the last decade where I've thought a scientist was wrong (or was failing to properly address a point I was making) and it turned out I was right. Every other time (including recently with Dr Judah Cohen, by email) it turned out I'd made a very basic beginners mistake.

So Occams razor suggests to me that the situation where amateurs think they've disproven the greenhouse effect (which has stood in its basics since Callendar and Arrenhuis), is most likely a case of '3'.

As I am busy, and would rather use the free time I have to study sea ice, I apply the following shortcut when dealing with such people: If you think you have disproven the greenhouse effect then publish your argument in the primary peer reviewed literature. Once the majority of physicists start to recognise that you are correct I will revisit the issue, until that point I am too busy to go over ground I have gone over in the past.

Once someone does disprove the greenhouse effect I will congratulate them as they receive their Nobel Prize for Physics. Because overturning an entire field is the sort of thing people get Nobel prizes for.

In the rare cases where an amateur physicist working as a patent clerk overturns a field of physics, recognition does come when the arguments are strong enough.

Chris Reynolds said...

PS

"As I am busy, and would rather use the free time I have to study sea ice.."

I hasten to add, this isn't my only hobby/pastime. I'm going to leave off films from here, but keep getting tempted to do some radio posts.

Mack said...

Thank-you Chris,
I do urge you to read and take consideration of the Nasif Nahle posts at Jennifer Marohasy's site.I'm of the opinion that you may be reading an historical "document" there. (there have been a number of attempts to sabotage her blog)
Luckily the internet never forgets. :)
Kind regards and all the best for 2015
Mack.

Chris Reynolds said...

Mack,

Hope you had a good New Year.

Anyway Nasif Nahle says in the abstract " ...the idea that a cooler system can warm a warmer
system are unphysical concepts".

Wrong, in a very basic way that's been put forward years before and rubbished back then when I was a sceptic and saw through it.

Consider a hot and a cold body. Both radiate photons. The net flow is from hot to cold. However just as there is a flow of photons from the hot to the cold, there is also a flow from the cold to hot. The net being the difference between the two. To suggest otherwise is to suggest that the photons 'know' the temperature they are being radiated from and know to steer clear of warmer objects. They do not!

Now lets say we raise the temperature of the cold body. It starts to radiate a greater flow of photons, and can in theory actually warm the hot body, albeit marginally so.

Actually, now I think about it I'm prepared to be really harsh. Nasif Nahle is talking utter bollocks.

Mack said...

"Now lets say we raise the temperature of the cold body"
Sorry...hypothetical, speculative nonsense Chris. The sun doesn't arbitrarily choose one body to radiate energy to, in preference to another. You need to stay real Chris.

Chris Reynolds said...

I am staying real Mack.

CO2 traps infrared radiation and leads to an increase in the temperature of the effective radiating layer of the atmosphere (high in the troposphere). This then has a knock on effect on temperatures in lower layers of the atmosphere and at the surface, which are warmer.

So the question is, can a warming of colder atmosphere above cause a warming at the surface? The physics show that it can.

The idealised experiment would be two black body emitters as parallel plates, in a vacuum. Plate 1 is hotter than plate 2. But if we raise the temperature of plate 2 we raise the emission power of plate 2. This increases the flux to plate 1, and plate 1's temperature must then increase in order to reach a new equilibrium point where emission is equal to the energy input warming the plate and the flux of infra-red radiation hitting the plate.

If you think Nsaif Nahle is correct you need to explain what happens to the energy increase from an increase in incident infra red as the cold plate is warmed. Or explain how the warmer plate does not intercept more infra red as a result of the increased temperature of the cooler plate.

Chris Reynolds said...

Actually there's too many words there. The core of the argument is whether increased back radiation due to CO2 coming from the colder atmosphere high above the ground can warm the ground.

jai mitchell said...

Chris,

Thought you would be interested in this NASA satellite animation of JJA insolation absorption anomalies. Note that the impact is from ice melting from the leading edges inward.

https://www.youtube.com/watch?v=YqQNrb8TvJE

Another thing to consider WRT cumulative insolation absorption during earlier and earlier months.

Chris Reynolds said...

Thanks Jai,

I hadn't seen that video, but it is precisely this consideration that has motivated me to make an index to better account for earlier open water under stronger insolation. This doesn't look like it is going to be much better than late July open water. However in the decades to come I think it will prove more useful than late summer open water.

R Graf said...

Chris,

Thanks for posting your Excel sheet. I am looking for a sheet or formula that should accurately model the seasons based on latitude and heat capacity. For example, the ocean seasons are muted by their increased heat capacity. This effect carries onto shore by wind and air convection and diminishes toward the land's heat capacity effect in accordance to the air's residence time over land.

I am looking not the reinvent a wheel if you have seen this anywhere.

BTW, regarding the GHE, I'll have a go at a logical way to explain it. Unlike cooling of everyday items on Earth, the only way the Earth itself can cool is by radiation since convection does not work in the vacuum of space. The sole interface the Earth has for outgoing infrared radiation is near the top of the atmosphere since water vapor, CO2 and other GHG blocks it any lower down. At the top CO2 is the major player in significant IR bands and since there is no vapor (below its freezing point). Higher CO2 concentrations narrow the effective radiation altitude range (closer to the top) for those IR bands. The closer one gets the the top the cooler. The cooler it gets the lower energy's efficiency to radiate. That loss of efficiency can only be made up by warming at every level below until the new steady state is achieved. It turns out that for every doubling of concentration of CO2 the temperature will rise 1.1C to establish the new steady state. This is accepted by 98% of scientist. What has scientists split is whether small warming gets amplified by clouds and other feedbacks or perhaps even gets muted by them. The IPCC best guess based on models is that amplification is about 2.8X due to increased atmospheric water vapor. So far this is unproven and the temp record for the last 130 years hints it might be only 1.5-2X. But if we get a huge jump this year (like .4C) as top climatologist Kevin Trenberth has predicted (and is praying for) then the IPCC models will not need to be scrapped because we have only had .1C rise in the last 17 years.

Chris Reynolds said...

R Graf,

I've just been looking but it seems all my notes on this have been lost between clearout and back up, so all I have is what is in the spreadsheet on Google Docs. This is annoying because I had a further formula accounting for optical dispersion varying with angle of the sun, the equation would have tied in with the term 'delta' in the above equations and would have reduced insolation at the polar regions.

Anyway at least you know that there is such an equation out there, not that this is much help. :(

Modelling the seasons is complex. In a maritime climate (like here in the UK) warming is tempered by the thermal mass of the oceans meaning that peak warmth is most typically in late summer, rather than near the maximum of insolation in early summer. In a continental climate the rise into summer warmth is more abrupt. But of course in winter the ocean tends to reduce the degree of winter cold in a maritime climate, whereas continental winters are colder. But one could go into greater complexity considering gyres in the oceans and impacts on heat transport...

The only way a simple numeric model could be done is by making some possibly rather uncomfortable simplifying assumptions. One could model a maritime climate over a single layer ocean, say 50m depth, and see the tempering of winter cold and the delayed warming into summer. What simplifying assumptions are reasonable would depend on the use to which you intend to put the model.

As for the GHE, my understanding is that the feedback on the basic radiative driven warming is positive, and that paleo data supports the Charney range (eg. Annan and Hargreaves). But I am having a break from sea ice until April and am too busy to get into a discussion on that: My day job (I'm an engineer) and some electronics I am working on in my spare time.

Anonymous said...

Hi Chris,
Comment on the average of your table - I think you may have forgotten to area-average? If you want the avg. for the globe, you need to take in to account the fact that 1 degree of latitude at the equator covers a much larger area than 1 degree near the pole. So you need to weight your average by the cosine of the latitude - this should pull your value of 295Wm-2 closer to 340 Wm-2 (=S_o/4. The total energy incident of the Earth's surface is s_0*area of Earth's disc. The average of this over the globe is s_0*area of Earth's disc/surface area of Earth. If you model Earth as a sphere, then this is s_0/4. In reality there is some small deviation as the earth is a slightly prolate spheroid).

Chris Reynolds said...

Hi Anon,

Perhaps, it's been nearly a year since I last looked at this and I'm too busy to get back into it. The table is just points at the vertices of the stated lat/lon, not a grid box average.

Anonymous said...

The maximum solar insulation is different by times of year due to eccentricity of the earth.

Chris Reynolds said...

Lulu,

Feel free to post a link to your improved model incorporating that and any other factors you want to incorporate.

Ed said...

Hey all,

I know this is years on from when this post was orginially made but to anyone like me who despaired when the orignial link to the excel spreadsheet didn't work I hope this alternative link helps.

https://drive.google.com/file/d/1CdmgrFNu9on2Otyvem7iGBEJy3EQSvVv/view?usp=sharing

This link will take you to a spreadsheet I made that calculates solar insolation on the factors given in this blogpost. The only inputs it requires is latitude and date.

Hope I have helped poor soul from struggling like I did.

(If you need help comment under the link and I'll try to respond)

Good luck!

Chris Reynolds said...

Thanks Ed,
Very much appreciated. I will try to find the original spreadsheet and link to it on my old Google Drive.
I have totally moved on from Sea Ice so everything has got rather tatty and ragged on here.
Chris

Anonymous said...

Hi Chris Reynolds,

Have you yet found your original spreadsheet?

This could be incredibly useful for me. No worries if not.

Thanks,
Richard

Chris Reynolds said...

Hi Richard,

I shall try to find it this coming week.

Chris Reynolds said...

Hi Richard,

I have updated the link in the blog post, and here it is: https://docs.google.com/spreadsheets/d/1-BSwCv8MEetnz_Ha2p9Ezckqrxk6TQXZ/edit?usp=sharing&ouid=116980323456999735226&rtpof=true&sd=true

Let me know if there's a problem with it.