A Study of Hadley-CRU weather station data

Do systematic effects caused by the geographic location of weather stations effect global temperature anomalies ?

I have been studying the station temperature data used in the Hadley Cru global temperature analysis. These data consist of all the land based stations dating back to about 1700 containing monthly averaged temperatures for each year of data. These are the data used to produce CRUTEM3 and HadCruT3. Hadcrut3 includes sea surface temperatures. I have been using and extending the PERL analysis programs kindly provided by Hadley – available here. As a by-product, I also developed a  geographic station browser allowing you to view and plot all individual station – for further details see also here. A comparison of the calculated anomalies based on the >3000 stations and the usual HadCRUT3 is shown below.

Fig 1: The Red points are the temperature anomalies used by IPCC - Hadcrut3V. The black points are the calculation results of station-gridder.pl | make_globl_averages.pl calculated from all the station data (as described below).

There are small differences which confirm that the main component of the observed temperature rise is based on the land data. However it is important to understand exactly how these data points are derived.

Algorithm:
1. Average monthly temperatures are calculated at each station based on the following procedure: record the minimum and maximum temperature for each day and then calculate the average of the minimum and maximum. Calculate the averages for the month from these daily data. These are then recorded for each station file along with metadata [lat,lon] station name etc.
2. So-called monthly “normals” are calculated for each station by averaging each individual monthly temperature over the years 1961 to 1990. Standard deviations are also calculated. The normals are then assumed to represent a “standard annual variation”.
3. Anomalies are defined for each station by subtracting the monthly values for a particular month from these normal values. Stations without normals for 1961-1990 or where any anomaly is > 5 standard deviations are excluded.
4. The world is divided into a 5×5 degree grid of 1592 points. For each month the grid is populated by averaging the anomalies of any station present within each grid point. Most grid points are actually empty – especially for those early years. Furthermore the distribution of grid points with latitude is highly asymmetric with over 80 percent of all stations outside the tropics.
5. The monthly grid time series is then converted to an annual series by averaging the grid points over each 12 month period. The result is a grid series (36,72,160) ie. 160 years of data.
6. Finally the yearly global temperature anomalies are calculated by taking an area weighted average of all the populated grid points in each year. The formula for this is $Weight = cos( $Lat * PI/ 180 ) where $Lat is the value in degrees of the midle of each grid point. All empty grid points are excluded from this average.

Quality of the data

Fig 2: Distribution of stations: Left are stations with data dating back to 1850: Right are stations dating back to at least 1930

From 1850 -1860 less than 5% of grid points contain data (figure 3). This rises to 20% by 1940 and peaks at 30% from 1960 to 1990, before falling again to  23% currently. Figure 4 shows the latitude distribution which demonstrates that over 80% of stations lie outside the Tropics at high latitudes ( Europe, US, Russia, Australia etc.). This can also be seen visually in the map displays of the Flash application (Figure 2).

Fig3: The percentage of the 5x5 grid containing any station data

Fig4: Number of grid points at high latitudes > +-25 compared to grid points in the extended tropics < +- 25. Less than 20% of the sample are in the (extended) tropical region.

Observations:

1. There is poor coverage over the main tropical warm zone of the tropics +-25 degrees. The averages are biased to high latitudes with large summer to winter swings. This then also accentuates the temperature differences between the southern and northern hemisphere.
3. The stations are all on land and lack the sea surface temperature measurements. However the annual anomaly results are almost the same as those from Hadcrutem3VG which include sea surface temperature data. Therefore the land based temperature data dominate the temperature trends.

Experiments:

The first excercise I did was to look directly at the temperatures rather than at the anomalies. This also shows how an unbalanced averaging a high latitudes accentuates differences between North and South hemispheres and the annual variations. The monthly temperature data for the full period is shown in Figure 5.

Fig 5: Global Monthly temperatures. The global average = (SH+NH)/2 . The smoothed curves are 50 point running averages

You can see how initially the discrepancy between north and south diminishes as more tropical stations are added. They then separate again after about 1920 before narrowing again recently. Note also how it appears to be the reduction in the minimum temperatures (e.g. for January – North)   that seems to be the cause of the anomaly rise.

Could there still be systematic effects due to changes in  the distribution and proportions of stations over time ? Next I looked at the averaged temperatures for each of 3 regions A) -20 <lat <20 (Tropics) B) Lat >20 (Northern lats) C) LAT<-20 (Southern Lats). These results are shown below. The year 1863 actually had no measurements inside the tropics – hence the zero value.

Fig 6: Yearly averages for Tropics, Northern and Southern Lats.

There are a couple of observations here. Firstly note  how excluding the tropics changes the “global average temperature” (B+C)/2 by only about 1 degree.  Secondly, note how most of the  temperature rise since 1980 is concentrated in northern latitudes. I am assuming that the “experts” choose to work with anomalies rather than  absolute temperatures because of known systematic problems. Anomalies measure the deltas between monthly temperatures relative to a standard(normal) set. In effect we are subtracting two large numbers from each other and averaging the residues. It is assumed that if the Earth is warming overall by say 0.5 degrees , then the global average of all the deltas measured by each individual station will also rise by 0.5 degrees. We are no longer measuring the global temperature as such but rather changes of an evolving distribution of station measurements over several decades. The anomalies can still in principal be prone to systematic effects through over-sampling. A simple example of  how this could happen would be say if North America rose by 1 degree while simultaneously the Sahara fell by 1 degree-  maintaining no net increase in global temperatures. The averaging algorithm would then produce a net global increase in temperature because the US is over-sampled while the Sahara is under-sampled.

Normalisation Study

Next I looked at possible effects of the normalisation method used to extract temperature anomalies for each station. The normal procedure followed by HAD-CRU is to calculate monthly averages for each month and each station between 1961 and 1990. These are then used to calculated anomalies for each station, and average them in each grid point. I decided instead to use the actual temperatures at each grid point resulting from the average of these stations. Then in a second step I calculate the monthly normals for each month at each grid point by averaging all the available data. There is no particular reason  to take a fixed time period for the normals, since  anomalies are just deviations from the norm.  First I generated temperature grids from 1850 to 2010. I then used the grid monthly time series to derive normals per month for each grid point. Then we subtract the normal from  grid temperatures to derive anomaly grids. Finally the area weighted and annual averages are derived. How does this compare with the standard result ?

Fig7: The blue points are temperature anomalies from GRID point averages compared to the standard anomalies based on per station based averages from 1961-1990 (black points)

As you can see the time trend changes dramatically! The results show clearly that prior to about 1910 the temperature anomalies are highly dependent on the normalisation method. This is not surprising as the geographic coverage is  < 14%.  After 1920, however this analysis shows that we can have a high confidence level in the published values.

Conclusions

1. The coverage of stations is concentrated at relatively high(low) latitudes. There is far less coverage in tropical regions, which thereby tends to exaggerates seasonal global temperature changes.

2. Sampling problems over in  5×5 grid for early years must lead to systematic problems in early years because there are so  many  empty grid points.

3. The standard use of temperature anomalies per station over a fixed period 1961-1990 is  somewhat arbitrary. When anomalies are instead calculated against a longer reference period and use grid values rather than station values then significant diferences are seen prior to ~ 1920.

4.The  conclusion is that the land based data are reliable after ~1920 but that earlier data are subject to systematic errors. The absolute temperature data are also effected  when there is significant asymmetries in geographic sampling.

 

Posted in Climate Change, climate science, Science | Tagged , , | Leave a comment

Global Temperature trends from stations post 1930

The map shows the location of those stations after 1930. “Mouse over” gives the name of the station and click shows the graph. You can zoom in by dragging a rectangle over an area. Please wait a couple of seconds to load !

Get Adobe Flash player

These are the stations used to generate the HADCRU3V dataset, and are kindly made available by the Hadley centre of the UK Met Office. For each station I have generated a time series graph showing temperatures and “anomalies” for the full time span. The Flash graph is by DIY map . The PERL GD package was used to create over 5000 graphs. Please note that GD can distort time scales because missing gaps of data are not shown. However trends in data are respected.

NOTES All station data for all time spans can be viewed in one application HERE

Posted in Climate Change, climate science, Science | Tagged , | Leave a comment

Global Station Temperatures beginning 1860-1930

The map shows the location of those stations with records back to between 1860 and 1930 (See next post for the others). “Mouse over” gives the name of the station and click shows the graph. You can zoom in by dragging a rectangle over an area. Please wait a couple of seconds to load !

Get Adobe Flash player

These are the stations used to generate the HADCRU3V dataset, and are kindly made available by the Hadley centre of the UK Met Office. For each station I have generated a time series graph showing temperatures and “anomalies” for the full time span. The Flash graph is by DIY map . The PERL GD package was used to create over 5000 graphs. Please note that GD can distort time scales because missing gaps of data are not shown. However trends in data are respected.

NOTES All station data for all time spans can be viewed in one application HERE

Posted in Climate Change, climate science, Science | Tagged , | Leave a comment

Global Weather Stations earlier than 1860

This is a tool for accessing long term weather data sourced by Hadley Centre and CRU. The map shows the location of those stations with data going back earlier than 1860 (See next post for the others). “Mouse over” gives the name of the station and click shows the graph. You can zoom in by dragging a rectangle over an area.

Get Adobe Flash player

These are the stations used to generate the HADCRU3V dataset, and are kindly made available by the Hadley centre of the UK Met Office. For each station I have generated a time series graph showing temperatures and “anomalies” for the full time span. The Flash graph is by DIY map . The PERL GD package was used to create over 5000 graphs. Please note that GD can distort time scales because missing gaps of data are not shown. However trends in data are respected.

NOTES:

1. All station data for all time spans can be viewed in one application HERE

2.  At one single location there is to my eye  no evidence of a systematic rise in temperatures locally. For example mean Dublin temperatures have not changed throughout continuous temperature measurements from 1850 to today. The results are shown below.  So it looks like we will still have to wait for a very long time before “Hawaiian” style warmth reaches the Irish Riviera!

DUBLIN Temperature and Anomalies relative to a 1960-1990 average

Posted in Climate Change, climate science, Science | Tagged , | 4 Comments

Sahara Temperature Trends

I have been interested in analysing temperature trends measured by weather stations in the Sahara region. This is for two main  reasons. The first is that with so little humidity in the atmosphere any effects caused just by increases in CO2 concentrations should be more evident. Water vapor is the Earth’s dominant greenhouse gas and any changes in  humidity either naturally  or perhaps linked to climate feedbacks  are difficult to untangle from CO2 increases. Consequently arid places with little changes in humidity should help to reduce this uncertainty. The second reason is that the Sahara has not seen such massive human development as seen in other places. I have downloaded the HADCRU station data released in  July 2011 from the Met Office available here.  I have also been using and modifying the  PERL files kindly provided which are used to calculated the gridded data for global temperature “anomalies”.  The details of this are taken from the HadCru site below.

“Stations on land are at different elevations, and different countries estimate average monthly temperatures using different methods and formulae. To avoid biases that could result from these problems, monthly average temperatures are reduced to anomalies from the period with best coverage (1961-90). For stations to be used, an estimate of the base period average must be calculated. Because many stations do not have complete records for the 1961-90 period several methods have been developed to estimate 1961-90 averages from neighbouring records or using other sources of data. “

To identify relevant Saharan stations for the study, I selected those stations whose latitude and longitude lie between LAT[15,28] and LON[-12,30]. When doing this I noticed that the files actually use Lon positive heading west which is oposite to normal so I had to correct for this. Each station was also required to have data covering at least the period 1960 – 2000. There were 18 stations which satisfied these criteria as follows:

IN SALAH      Algeria
TAMANRASSET   Algeria
BILMA         Niger
AGADEZ        Niger
TESSALIT      Mali
KIDAL         Mali
TOMBOUCTOU    Mali
GAO           Mali
NIORO DU SAHEL Mali
NARA          Mali
HOMBORI       Mali
MENAKA        Mali
BIR MOGHREIN   Mauritania
TIDJIKJA      Mauritania
SEBHA         Libya
KUFRA         Libya
DAKHLA        Egypt
FAYALARGEAU   Chad

For each of these stations the data provide the average monthly temperature over time and the calculated “normals” and standard deviations. Using these temperature values and after subtracting the “normals”, the anomalies for each of the 18 stations are then calculated. An example “snippet” from the station file for FAYALARGEAU, Chad is shown below

Number= 647530
Name= FAYALARGEAU
Country= CHAD
Lat=   18.0
Long=  -19.1
Height= -999
Start year= 1946
End year= 2007
First Good year= 1946
Source ID= 10
Source file= Jones+Anders
Jones data to= 1977
Normals source= Data
Normals source start year= 1961
Normals source end year= 1977
Normals=  20.2  22.5  26.1  30.6  33.0  34.0  33.1  32.7  32.6  29.4  24.6  21.3
Standard deviations source= Data
Standard deviations source start year= 1946
Standard deviations source end year= 1977
Standard deviations=   1.6   1.8   1.6   1.2   1.2   0.9   0.9   1.0   1.2   1.3   1.6   1.6
Obs:
1946  21.9  21.2  26.4  31.8  35.3  35.4  34.3  33.7  34.1  31.2  27.2  22.3
1947  20.9  26.7  26.8  31.2  34.4  36.0  35.0  34.5  34.0  30.9  25.2  23.3
1948  21.4  22.8  23.7 -99.0  34.9  35.4  34.2  33.9 -99.0  28.6  23.7  19.2
1949  20.2  19.7  26.3  28.9  34.0  33.5  32.5 -99.0  31.3  28.8  24.7  19.0
.......

To summarise: The calculation of the grid of global average temperature anomalies used by Hadley-Cru is based on a per station “normal” set of monthly data values and associated standard deviations for 1961-1990.  These normal  values for each station are the monthly averages for that particular station over the period from 1961 to 1990. For the Saharan study the 18 station anomalies are all plotted together in Figure 1. There appears to be a very slight rise at recent times above the red zero line.

Fig 1: Plot of all 18 Station anomalies. The red line is the zero offset to the monthly normal values.

I next take the temperature values for each of the 18 stations and average them all together to get a single time series. This is shown in Figure 2.  together with a 50 point rolling averaged smoothing term. There is little evident change in mean temperatures. Figure 3 shows the averaged  anomalies for all 18 stations for the full time period, together with a least squares linear fit, and a clear rise in net anomaly now becomes evident. The trends between the averaged temperatures and the averaged anomalies  are significantly different. The averaged anomalies as used by the HadCru  Gridding algorithm shows an almost linear slow rise of about 1 degree C.  in temperature anomaly from 1950 to 2011. However there is no real sign of this at all in the station averaged temperatures over the same time period. Instead it is just the temperature extremes – both  positive and negative which  seem to increase slightly.

Fig 2: The Average monthly temperatures for the 18 Stations. Shown too is a 50 point smoothed average which shows little change in "average" temperature although the monthly extremes seem to increase.

Fig 3: Average of the calculated anomalies from the 18 Stations using provided monthly normals. Shown is a linear least squares fit through the data

Could there be a systematic effect in the analysis which accentuates any small increase when averaging together all the temperature anomalies ? To investigate this I calculated a new set of anomalies based on all 18 stations. The normals were calculated by averaging temperatures over all 18 stations for each month over the period 1961-1990. The “net” anomalies were then calculated by subtracting these “averaged” normals from the averaged temperatures over the full period. The result is shown in Figure 4.

Fig 4: Anomalies for the average temperatures for 18 Saharan Stations. The red line is a least squares fit through the data.

The anomaly trend is now essentially flat. What can be the difference between using the two different normalisation values ? It seems that averaging the anomalies gives a slightly different result to the anomalies of the averages (figure 5). However the conclusion is crucial as to whether temperatures have risen overall in the Sahara over the last 60 years. In the first case the conclusion is that average temperatures in the Sahara have increased since 1950 by just over 1 degree.C while  CO2 concentrations increased by about 80 ppm – in line with IPCC predictions.  In the second  case the conclusion is that there has been no significant change in temperature at all during the last 60 years. The raw monthly average temperature measurements also support little or zero net change. I have always worried about the use of “anomalies” instead of temperatures because it assumes there is one “normal” period (1961-1990) to which all other temperatures should be referenced. Could this assumption itself introduce a bias ?

Fig 5: Detailed Comparison between the two anomalies. Both are normalised to the period 1961-1990.

Posted in Africa, Climate Change | Tagged , , | Leave a comment

Observation of a Lunar Signal in Solar Irradiation Data

Late last Friday evening as I wasted a couple of hours trying in vain to buy 2 tickets from the London  Olympics website, I filled my time looking further at the SORCE TIM solar irradiation data [1]. The olympics website turned out to be broken but all was not in vain because instead, I had a brainwave about what causes the apparent regular small oscillations in the data. It is the Moon of course ! In reality the Moon and  the Earth orbit each other around  their common centre of mass. This is located about 4000 km outside the center of the Earth so as a result the Sun-Earth distance also  changes in phase to the lunar month. The exact monthly  variation in distance (D) will depend on  orbital parameters but could be as much as 8000km.   The data clearly show a monthly variation of about 0.2 watts/m2  (figure 1) during 2008-2009 which turns out to be  in  phase with the lunar cycle. There is a negligible probability that this occurs by pure chance – see below.

Fig 1: The left axis shows the TIM data and the theoretical data for Earth-Sun orbit. The difference between the two curves is plotted as the black signal scaled to the right y-axis. The red curve is a calculation of the lunar phases based on a siderial month of 29.53 days using an arbitrary scale for comparison. At full moon the Earth is approximately 4000km nearer the sun.

As far as I am aware this is the first time this effect has been seen as  I cannot  find any other reference to  in the literature. Apologies if I am mistaken.

My original aim in studying the TIM data was to try to isolate a possible (annual)  signal as the Earth passes through plane of the interplanetary dust cloud. Could ice ages be really  caused by  an oscillation in dust density lying between the Earth and the Sun in phase with the eccentricity of the earth’s orbit ? But that is another story…

References

[1] http://lasp.colorado.edu/sorce/index.htm

Posted in Climate Change, Ice Ages | Leave a comment

Evidence of variations in solar irradiance with Earth’s orbit

Presented here is evidence of variations in solar irradiation as the Earth follows it’s elliptical orbit around the sun. It is proposed that these variations are due to changes in the density of dust clouds lying between the Earth and the Sun as the orbit moves from perihelion to aphelion. These variations become evident when compared to a calculation of the expected irradiation based on the Earth’s orbital parameters and a measured average solar constant of 1361 watts/m2. The change in irradiation between perihelion and aphelion is about 2 watt/m2 for the current eccentricity of 0.0167. This supports the hypothesis that the primary cause of recent Ice Ages has been  increases in interplanetary dust for lower eccentricity of the Earth (“resonant” clouds).  If the difference between an expected 1368 watts/m2 average solar irradiation and the measured value of 1361 watts/m2 is assumed to be all due to absorption and scattering of dust lying between Earth and the Sun then current absorption varies between a maximum of 8 watts/m2 at perihelion(closest approach) and a miniumum of 6 watts/m2 at aphelion (largest distance). During long glacial periods average global temperatures fall to about 2 degrees below current values which infer a further drop in solar irradiation of about 7 watts/m2. Therefore in low eccentric orbits the solar constant would then fall to  ~1354 watts/m2. The hypothesis is that for the last million years interplanetary dust resonates to the same same gravitational cycles as Earth in the ecliptic plane leading to higher densities at lower eccentricities. This is described in the previous post

The data is taken from  the SORCE Total Irradiance Monitor (TIM) [1] whose long-term point to point uncertainty is estimated to be < 0.014 W/m^2/yr (10 ppm/yr) , and has an absolute uncertainty of 0.48 watts/m2.  The calculation of the solar intensity  uses a solar constant of 1361 watts/m2, an eccentricity of 0.0167 and tracks the orbit for each day applying conservation of angular momentum and a year of length 365.25 days. The comparisons are shown in figure 1. Also plotted is the difference beween calculated and observed irradiance (first y-axis).  A clear regular discrepancy is observed with some high frequency substructure evident as well.

Fig 1: Comparison of Measured Solar Irradiance from SORCE TIM Total Solar Irradiance instrument with calculated

Fig 1: Comparison of Measured Solar Irradiance from SORCE TIM Total Solar Irradiance instrument with calculated

The amplitude of the effect is about 2 watts/m2. At perihelion irradiation is attenuation greater than expected while at aphelion irradiation is greater than expected. This can be understood assuming a constant average attenuation of X watts/m2 with X+1 and X-1 representing net changes in density. Stefan Boltzman’s law predicts S0=1368 watts/m2 which all else being equal implies X=7 watts/m2.

A detail of a single orbit is shown in figure 2.

Fig 2: Detail of one orbit showing higher irradiance at aphelion and lower irradiance at perihelion as compared to vacuum values. Note the slight offset of the signal.

Fig 2: Detail of one orbit showing regular small variations compared to vacuum values. Note the slight offset of the signal.

If we assume that the 2 watts/m2 discrepancy is due only to absorption by intermediate dust, then the density between the Earth and the Sun varies by about 30%. For comparison the total estimated radiative forcing caused by all CO2 emissions since the industrial revolution is also about 2 watts/m2.

References
[1] http://lasp.colorado.edu/sorce/index.htm

Posted in Climate Change, Ice Ages | Tagged , | Leave a comment

Part 2: The real cause of Ice Ages ? – Resonant dust clouds ?

This post proposes that regular variations in interplanetary dust between the Earth and the Sun is the primary cause of  recent ice ages whose oscillations have been parameterised here. The Milankowitz theory of insolation is able to well explain the continuous 41,000y climate oscillation but it cannot explain why a 100,000y cycle became prevalent during the last million years, even though it is in phase with orbital eccentricity of the Earth’s orbit around the sun. A different origin for the emergence of the 100,000y cycle of major ice ages is proposed. This is based on the hypothesis that a large comet broke up close to Earth’s orbit about 3.3 million years ago, leading to an increase in inter-planetary dust whose density eventually phase locked to variations in Earth’s eccentricity. This process is able to explain the cooling trend starting 3 million years ago eventually leading  to major eccentricity driven glaciations with persistent modulation in insolation from the obliqueness cycle. Evidence in support of regular cycles in inter-planetary dust density have been found in ocean sediments, and there is some evidence of a comet fragment impacting with Earth at 3.3 million years ago.

Insolation

Milankowitch cycles inducing insolation changes in the polar regions is the generally accepted primary cause of Ice Ages. However the details for insolation changes due to the ellipticity of the Earth’s orbit don’t really work out. Although the insolation varies by 6%(today) up to  20% between perihelion(currently southern hemisphere summer) and aphelion(currently northern hemisphere summer), the net annual change for the Earth is insignificant. The poles themselves are unique places because they essentially experience just 1 day and 1 night per year. For the North Pole the sun rises at the spring equinox and sets again at the autumn equinox. At mid summer’s day the 24 hour insolation reaches average values higher than at the equator. The winter night then follows resulting in  zero insolation for the next 6 months. This also explains the rapid  melting and freezing cycles observed in polar regions.

Data used for insolation calculations at the poles are the three Milankowitz cycles of Orbital Eccentricity, Obliqueness and Precession originally calculated by [1] are taken from the Nasa site[2]. The insolation is then calculated by tracking the orbit through the year on an hourly basis. Sample results are shown below.

Variation in insolation at the pole for 3 different (North=Perihelion, South=Apehelio, Equinox=mid point), and obliqueness (Phi)

Variation in insolation at the pole for 3 different alignments (North=Perihelion, South=Aphelion, Equinox=mid point), and 3 different obliqueness (Phi)

The effects on maximum insolation, which  effects summer melting rates  of ice at the poles, is evident for both obliqueness and precession. However the effect of changes in the orbital ellipticity is more subtle. An increase in insolation for one pole close to perihelion is offset by a diminished insolation for the other pole. Later in the precession period  when the poles are positioned at the equivalent of todays equinoxes, the polar insolation  is similar to that of a circular orbit. Overall the net annual effect of eccentricity of the Earth’s insolation is negligible.  For the last 900,000 years however,  Ice Ages show a clear correlation with ellipticity at 100,000 years and this remains unexplained. The basic problem can be seen better by calculating long term insolation trends. Figure 2 shows the maximum and total insolation for the poles over the last 600,000 years.

Maximum and total solar insolation at the pole during last 600,000 years

The main effect of ellipticity is to modulate the effect of the precession terms at 23,000 years. Indeed if the orbit was circular then the precession of the equinoxes would have no effect whatsoever. However the direct effect on insolation of the increase in ellipticity of the Earth’s orbit averaged over a year is negligible. This is the paradox as to why regular glaciations and warm interglacials occur every 100,000 years in phase with the eccentricity.

The clearest insolation effect on climate is the obliqueness signal which when averaged out over the combined ellipticity and precession yields a regular variation of around 25 watts/m2 in maximum flux. The delta 18O data show a clear signal throughout the 5 million year period in phase with obliquity. The precession signal is directly dependent on eccentricity and is less evident in the data. However there is no explanation for the primary 100,000 year variation in phase with ellipticity. An excellent review of this problem has been given by  Maya Elkibbin & JoseA. Rial [3] which highlights the questions that any complete theory of the evolution of Ice Ages must answer namely:

1. The observed gradual decrease in temperatures starting just  over 3 million years ago.

2. The 43,000 year continuous signal dominated until about 1 million years ago when the onset of  70-100,00000 glacial cycles begin. Why did this happen ?

3. What is the cause of the recent 100,000y glacial cycles, since Insolation changes caused by ellipticity are too small to be an explanation ?

4. Why did the frequency change from about 70,000 years between 900-700,000 ybp to the current cycle of 100,000 years ?

5. Why is the larger 400,000 year eccentricity signal absent from recent data ?

Variations in insolation can explain the 43,000y obliquity and 23,000 y precession signals observed in the data. However another explanation is needed to explain the points made above. A possible solution is proposed as described below

Interplanetary Dust.

Muller et al. have proposed an alternative explanation for the 100,000y cycle. They suggested that the driver is the change in inclination of the Earth’s orbit. This follows an 80,000 year cycle relative to the current ecliptic plane, but when projected to the invariate plane defined as the centre of mass plane of the solar system this is extended to 100,000 years. The physical process envisaged is that the interplanetary dust cloud is centerd on the invariate plane and the Earth passes through this during its orbit. Evidence from ocean sediments show a 100,000 year cycle of H3 deposits associated with cosmic dust [3]. In effect, changes to orbital parameters of eccentricity, inclination and longitude of ascending node are all caused by the mutual gravitation of solar system planets as they orbit the Sun. Similarly all bodies in the solar system will follow such cycles including dust clouds. Let’s investigate correlations of inclination and glaciations. Figure 4 shows the data as given by [5] compared to the benthic foram stack.

Comparison of data with inclination and eccentricity signals

Comparison of data with inclination and eccentricity signals

The hypothesis that the Earth passes through a static interplanetary dust cloud aligned with the invariate plane  attenuating sunlight is attractive. Furthermore there is direct evidence in ocean sediments of a regular 100,000 cycle of He3 deposits [4].  Unfortunately the data do not seem to support inclination as being responsible, despite the similarity in timespans. One would expect colder temperatures at low alignment where presumably the Earth crosses the dust plane. The correlation with eccentricity is better with high eccentricity seemingly algning with warmer interglacial periods. So what is really happening and can we explain how the climate changed to regular glaciations ?

First we try and estimate if there is today an observable effect of interplanetary dust attenuating the  solar constant as measured accurately by satellite.

The surface temperature of the sun is 5778 K and the mean distance of Earth to the Sun is 1.496×108 km. Applying Stefan Boltzman’s law we can calculate the radiation flux (solar constant) out side the atmosphere. Radius of the Sun = 6.96*10^4km

S0 =  (sigmaT^4)*4piR^2/4pi*D^2   putting the numbers in gives:

S0 = 1368 watts/m2

The measured average solar constant over the orbit is 1361 watts/m2 [6]. This implies therefore that currently during an interglacial maximim on average 7 watts/m2 is absorbed (or scattered) by interplanetary dust currently lying between the earth and the Sun. In addition accretion of dust by the atmosphere which is estimated at 40000 tons per year will also slightly reduce sunlight incident on earth. Therefore it is proposed that variations in the dust cloud aligned with the Earth – Sun ecliptic plane is the origin of recent Ice Ages. Until 1 million years ago, colder glacial periods were driven by variations in obliquness induced insolation. However this changed presumably when the dust cloud settled into harmony with the Earth’s orbit.  For the last million years “Milankowitch” variations in the eccentricity of the orbit have also applied  to the dust cloud inducing regular variations in density distributions thus causing variations in the solar constant. A hypothesis that can explain the slow cooling of the Earth, the onset of ice ages and their transition  from 43,000 y to 100,000y dominant cycles  is as follows.

Some 3.3 million years ago a large comet breaks up near the Earth creating a large irregular dust cloud lying between the earth and Sun. There is some direct evidence of a collision  of a significant cometry fragment with the Earth at this time from Argentina [5]. Initially this irregular dust cloud orbited the sun and acted to gradually reduce temperatures on Earth as it condensed into a Earth-Sun ecliptic orbit. The cloud gradually spread out into a disk in synchronous orbit with the Earth around the Sun. Variations in the  orbital eccentricity of the Earth are caused by the regular combined gravitational forces  of the outer solar system planets. Exactly the same forces will apply to the dust cloud now in synchronous orbit. However,  the low mass of particles in the cloud accentuate the orbital effects especially in eccentricity. Thus we can envisage a current dust cloud whose density and shape is moulded by those same eccentricity forces as that  on the Earth. Initially the diffuse cloud caused global cooling for 2 million years after breaking up. Slowly as the dust cloud settled into the ecliptic plane it  became resonant with the orbital effects of the Earth. When the Earth is at a maximum of eccentricity in the 100,000y cycle the dust cloud is spread in an elliptical cloud with much of it lying outside the Earth Sun line of sight. As the Earth resorts back to low eccentricity and a more circular orbit so too does the dust cloud leading to an increase in dust density lying between the Earth and the Sun. This effect is shown schematically below.

Effect on interplanetary cloud density with orbital eccentricity of the Earth for a cloud resonant to the Earth-Sun orbit.

Schematic of interplanetary cloud density variation with orbital eccentricity of the Earth assuming a cloud resonant to the Earth-Sun orbit with a similar planetary gravitation forces.

If it is assumed that the net effect of the current dust density between Earth and the Sun is such as to reduce the solar constant by 7 watts/m2, as described above,  then a doubling of this density as eccentricity drops to a minimum would cause a global temperature drop of about 2K leading to a new  an Ice Age. This is based on a Stefan Boltzman respons of :  DS = 4sigmaT^3 DT

Predictions

This model of an interplanetary  dust cloud resonant with the Earth’s eccentricity cycle predicts that the following phenomona should be observable.

  1. A slow decrease in the solar constant as the density of dust increases with falling eccentricity.
  2. Evidence of more dust external to the earth’s orbit along the aphelion axis  compared to the aehlion axis.
  3. A higher density of dust between the Earth and the Sun at perihelion rather than at aphelion. This could be reflected in a slightly lower value of solar constant at perihelion than otherwise predicted. In addition the Zodiacal light should be more intense at perihelion than aphelion.

AGW

A doubling of CO2 levels in the atmosphere by the end of this century will cause a direct radiative forcing of 5.3Ln(2) watts/m2 or 3.67 watts/m2. IPCC models usually assume a further feedback mainly from water vapor of 2 watts/m2/degC. A reasonable assumption would be that AGW would induce a warming term of about 5-6 watts/m2 after which CO2 levels are expected to  stabilize. The net effect of all human induced global warming could be to increase radiative forcing by something like between 4 and 8 watts/m2.

The eccentricity of the Earth’s orbit is currently decreasing with cooling beginning in 2000 years time leading eventually to another intense Ice age after  20,000 years and lasting a further 60,000 years. A new Ice Age would be disastrous for human and animal life at large latitudes.  If this theory of a variable interplanetery dust cloud is correct, then  AGW will likely  delay the early  onset of the next  ice age. To fully offset the cooling  effects of a thicker dust cloud, it would be necessary to keep CO2 levels above about 700 ppm for another 70,000 years.

Summary

A new possible cause of  Ice Ages is proposed based on an increase in interplanetary dust for lower eccentricities of the Earth’s orbit for the last million years leading to large cooling of the Earth.  This is proposed as the cause of the switch from  glaciations driven by the 43,000 year variation of insolation due to obliquity to a 100,000 eccentricity driven cycle.  The next cooling cycle should begin in 2000 years time reaching a glacial minium in about 10,000 years time lasting 70,000 years. At the same time anthropogenic global warming may delay the cooling trend if induced temperature rises of 2-3 degrees occur and defer the next glaciation if CO2 levels were to stabilize at ~700 ppm for 80,000 years.

References

[1] T.R. Quinn, S. Tremaine, & M.A. Duncan A three million year integration of the Earth’s orbit, Astronomical Journal 101, 2287Ð2305, 1991

[2] http://aom.giss.nasa.gov/srorbpar.html Andre L. Berger, 1978. Long Term Variations of Daily Insolation and Quaternary Climatic Changes. Journal of the Atmospheric Sciences, volume 35(12), 2362-2367

[3] Maya Elkibbi), Jose ? A. Rial, An outsider’s review of the astronomical theory of the climate: is the eccentricity-driven insolation the main driver of the ice ages? Earth-Science Reviews 56 Ž2001. 161–177

[3] Higgins et al. Sediment focusing creates 100-ka cycles in interplanetary dust accumulation on the Ontong Java Plateau, Earth and Planetary Science Letters 203 (2002) 383-397

[4] P. H. Schultz, M. Zarate, W. Hames, C. Camilión and J. King, A 3.3-Ma Impact in Argentina and Possible Consequences, Science 11 December 1998: Vol. 282 no. 5396 pp. 2061-2063

[5] Surface Temperature of the Sun http://en.wikipedia.org/wiki/Sun

[6] Kopp, G. and J. L. Lean, A new, lower value of total solar irradiance: Evidence and climate significance, Geophys. Res. Lett., 38

Posted in Climate Change, Ice Ages | Tagged , , | 1 Comment

What causes interglacials? Part 1.

The natural climate of Earth today is really an Ice Age with temperatures 2-4 degrees colder than today. For the last 1 million years glacial periods lasting about 80,000 years dominate interspersed with regular warmer interglacial periods lasting 10-20,000 years. The rapid expansion of human civilization has all occured during the most recent warm period which has already lasted about 10,000 years. In the previous post it was shown that each interglacial coincides with a maximum of the orbital eccentricity. Furthermore regular smaller temperature rises also occur with larger obiquity. The last interglacial occured with a coincidence of both effects which may have caused the rapid rise in temperature from the last glacial maximum.

The strange thing about both these  observations is that each effect acts to accentuate differences of insolation with season. Currently the Earth’s ellipticity is 0.0167 and summer in the Southern Hemisphere coincides with the closest distance to the sun. This also means that winters in the Northern hemisphere are less severe now than they were at the last glacial maximum when 13,000 years ago due to the precession of the equinoxes.

Ellipticity of Earth's orbit round the sun

Ellipticity of Earth

We can calculate the change in average insolation during the year with ellipticity assuming an average value of 342 watts/m2 for a circular orbit. The result is shown in the next figure.

Change in solar radiation during the year for different ellipticity

Change in solar radiation during the year for different ellipticity

Conservation of angular momentum ensures that the average solar radiation remains 342 watts/m2, but this also shortens slightly the seasons can be seen in the variation of cross-over for different ellipticities. The value 0.058 is the largest ellipticity the earth reaches during a 400,000 year cycle. The current value is shown in red and shows that the southern hemisphere summers receive up to 25 watts/m2 more radiation than Northern summers and the inverse for winters.  What does the actual data show ? We can see this by taking the Hadley Cru global average temperature data since 1850. Instead of looking at global warming we just look at the comparison between the monthly temperature anomaly versus the yearly average temperature anomaly. Then look at he average divergence over the 160  year period. We get the following curve.

Temperature divergence for each month for Full Hadcrut3

Temperature divergence for each month for Full Hadcrut3

The result shows an effect of around +- 0.04 degreesC but in the wrong direction !  This shows that the primary effect on Earth is the diference in ocean coverage between North and South. Despite the southern hemisphere receiving more solar energy in summer, this is offset by the thermal inertia of the oceans. Greater temperaure extremes are seen in the Northern hemisphere, than in the southern hemisphere. For this reason Ice Ages and interglacials seem to be triggered by effects mainly in the Northern Hemisphere

It is clear that the ellipticity effect is somehow getting magnified in the Northern Hemisphere to bring each Ice ages to an end. In the tropics the variations in solar energy effectively gives rise to an extra orbital season. 17,000 years ago thye variation in insolation at the equator was 10% between equinoxes due to this effect. There is also a subtle effect on global average temperaures caused by the effective negative feedback from Stefan Boltzman. Differentiative Stefan Boltzmans relation we get

ds = 4*sigma*T^3 dt    so  a temperature change dt will be induced by a change in radiation ds at temperature T of  ds/(4*sigma*T^3). This can be calculated  relative to say a fixed T=283K for a circular radius averaged over the year. This results in a small warming effect for large eccentricities shown below.

Average Temperature versus ellipticity for different feedback F

Average Temperature versus ellipticity for different feedback F

When the Earth’s orbit is near circular then summer and winter are symmetric between hemispheres. The Earth’s natural temperature would appear to be some 2-5 degrees cooler than today and in this state glaciers build up in Northern latitudes. As the orbital parameters change then increases in solar radiation for northern summers begins to melt the ice and lower albedo. This then seems to cause a run away effect which takes the earth to a new state of warmer interglacial period. The last ice age came to an abrupt end when both the ellipticiy and obliquity were at a maximum AND northern summers coincided with maximum insolation. After a warm period so the cycle repeats as the orbit becomes more circular again and the tilt of the Earth lessens and ice begins to build up again.

A fit to 800,000 years of glaciation cycles based on changes in ellipticity and obliquity was described in the last post. The overall agreement is remarkably good and shows that naturally the Earth will begin cooling to a new Ice age within 2000 years. Another Ice age would be dissasterous for mankind and it is difficult to see how the Earth could support the current population. Anthropogenic Global Warming (AGW) resulting from enhanced greenhouse gasses will likely cause between 1 and 3 degrees of warming within the next 100 years after which the temperature should level off. All this depends on climate sensitivity (feedbacks from water). When seen in the context of a million years of glacial cycles such warming could turn out to have one beneficial effect, in the sense that it could delay the inevitable astronomic cycle leading to the next Ice Age. However, it is naive to imagine that Man can overide an inevitable astronomic origin for glacial cycles.

Posted in Climate Change, Physics | Tagged | Leave a comment

Phenomenology of Ice Ages

The exact cause and evolution of ice ages still remain a mystery [1]. Why did the Earth shift into an unstable climate with global oscillations of 4-5 degrees every 100,000 years? What is the cause of the large scale glaciations with shorter warm inter-glacials? Proposed Milankowitz cycles cannot properly explain why their relatively small radiative effects apparently have been amplified during the last million years or why long term temperatures have been falling for over 20 million years. The work presented here does not attempt to explain the origin of  the ice ages, but instead it tries to  parameterise the observed temperature dependencies to derive as much quantitative information as possible. This approach applied in other branches of physics is called phenomenology[2]. The amplitudes of Milankowitz harmonic cycles are derived from making fits to data.  These fits are then used to predict the occurence of the next Ice Age with some confidence.

The “LRO4″ data [3] is a remarkable set of global proxy temperature data spanning a 5.3 million year period throughout the pleistocene Ice Ages.  It consists of a  stack of 57 globally distributed benthic d18O records. Unlike the Vostok Ice core samples which measure temperature at a single site in Antarctica, LRO4 represent the global temperature  and should therefore be more representative of climate changes. The data are based on samples of deep sea sediments consisting of calcium carbonate from plankton. The oxygen content (CaCO3) consists of a small amount of O18 mixed with O16. The ratio of O18 to O16 increases slightly with cooler surface temperatures. This is because O18 is slightly heavier and evaporates less readily than O16. Delta O18 is the ratio of  O18/O16 compared to a standard Calcite value. LRO4 represents a global measure of this ratio from ocean sediments and can be directly related to prevailing temperatures at the time. It provides a unique overview of the onset of Ice Ages (Pleistocene Era) starting 5.3 million years ago.

The objective of this post is to try and quantify the size of Milankowitz components in the LR04 data by fitting the data to a linear sum of Milankowitz harmonics plus any overall linear trends. The data as presented fall into 5 natural time periods based both on sampling rate and clear differences in signals. These are 1) 0-600 thousand ya 2)600-900 thousand ya 3)900-1500 thousand ya 4) 1500-3000 thousand ya and 5)3000-5300 thousand ya. For each period  a fourier transform was used to identify any harmonics present and then the data was fitted to a general  form with up to 3 harmonics and a linear offset with optional slope

deltaO18 = const + ax + b*sin(2pi/wavelength*x – c) + d*sin(2pi/wavelength2 – e)+ f*sin(2pi/wavelength3 - g)

Results:

1) For the last 600,000 year data there are clear oscillations present at 100,000 years and 41,000 years corresponding to the eccentricity changes of the Earth’s orbit and the change in obliquity. A very small signal for 23,000 years was seen but including it in the fit made little difference to Chi-squared. It seems to have a negligible effect on Ice Ages. The best fit result was:

delO18 = 4.0 +0.407*sin(0.0628*x-39.5) + 0.259*sin(0.153*x-45.3)

2) 600 – 900 kya

Already at 600 kya there is a change evident in the cycle of Ice Ages. The 41,000 year cycle is clear but the time delay between successive Ice Ages has begun to shorten from 100,000 years. In fact a fit using 100,000 years is not good.  A much better fit is found with a time period of 84,000 years. The reason for this change in the data is not clear. No 23,000 year signal is present. The best fit found was

delO18 = 4.2 +-0.36*sin(0.075*x-31.1) + 0.29*sin(0.153*x+5.27)

The combined fits compared to the data from both eras combined is shown below.

Combined fit back to 900,000 ya

Combined fit back to 900,000 ya

3) 900 – 1500 ya

Big changes begin to be evident in the data further back than 900,000 years. The data set from 900-1500 shows that the 100,000 year cycle has almost disappeared leaving just the 41,000 year signal. The best fit found was del18 = 3.8 + 0.1*sin(0.0628*x+19.7) + 0.248*sin(0.153*x+11.0). There is then a clear linear warming trend from 1500 ya until 3000ya where the only evident signal is that from the oblicity (41,000 year cycle). The best fit was a) from 900-1500

3.8+0.1*sin(0.0628*x+19.7)+0.25*sin(0.153*x+11)  - (small 100,000y cycle)

and b)

del18 = 4.0 – 0.00049*(x-1500) -0.17*sin(0.153*x-4.5) – 0.006*sin(0.273*x-29.7)

The third term is the 23,000 year precession term and the tiny amplitude demonstrates that it has a negligible effect. The combined fits from

Combined fits compared to data 900-1500 ya

Combined fits compared to data 900-3000 ya

Finally the fit to data from 3000 to 5300 kya begins to show a longer 400,000 year eccentricity oscillation as well as t 41,000 year oblicity. The best fit found was

delO18 = 3.25 -0.0002*(x-3000) -0.03*sin(0.015*(x-3000)-49.6) +0.07*sin(0.153*(x-3000)-31.3)) and is shown below.

Evidence of a 400,000 year elipticity oscillation 3000-53000 ya

Evidence of a 400,000 year elipticity oscillation 3000-53000 ya

The full 5.3 million years of global temperature change as reflected in the LR04 data with the fits described above is shown in figure 5.  This clearly shows a long term cooling of the climate leading  to the severe Ice Ages of the recent past. Three proposals have been given for this long term cooling trend. The first is that Antarctica has moved to fully cover the  South  Pole, isolating any heat flow from the oceans and lowering the Earth’s albedo [4]. The second proposal is that the Tibetan plateau has been slowly rising as India continues to crash into Asia[5]. The increasing ice cover over a large tropical region then changed the Earth’s energy balance by  decreasing albedo.  The third proposal is that the closing of the Panama Isthmus 3 million years ago altered the flow of heat in the oceans leading to the Gulf stream[6]. This brought more moist air to the polar regions resulting in more snow and ice.  The Earth’s climate appears currently to be in a highly unstable state, with sawtooth-like oscillations every 100,000 years flipping between long glaciations and shorter warm inter-glacial periods. The onset of a new Ice Age is a gradual process ending when temperatures reach a lower limit flipping  back rapidly to a warmer phase for 10-20,000 years. The fits demonstrate that the driver for this instability are the two cycles – eccentricity and obliquity. However the details of the correlation are surprising.

For the last 600,000 years Ice Ages have followed a 100,000 year cycle coinciding with the Milankowitz cycle. Every inter-glacial period has been in phase with a maximum in the eccentricity. This corresponds to increases in elipticity of the orbit, which seems counter-intuitive since it accentuates the variation in solar energy throughout the year. One fact though is that a smaller perihelion increases the maximum solar insolation. Unlike Ice core data the LR04 data are global and therefore must represent changes in total energy balance.  From 600-900 million years ago the duration of Ice Ages shorten slightly to 84,000 years and further back disappear completely. It is as if the long term cooling trend needed to reach a threshold for the orbital eccentricity oscillation to kick in.

The second clear observation is  that the 41,000 change in the angle of rotation of the Earth to the orbital plain (oblicity) has had a major effect on climate for at least the last 5 million years.  Surprisingly, the larger the tilt the warmer the climate, as demonstrated in figure 4. A larger obliquity means there are larger extremes in summer and winter temperatures. Both the Arctic Circle and the tropics increase in extent.  I suspect it must be the second effect which wins out, and when this is coincident with  large eccentricity with an annual closer distance to the sun – Ice Ages end !

Correlation of larger obliquity and warmer temperatures

Correlation of larger obliquity and warmer temperatures. Note that lower deltaO18 mean warmer global temperatures. The fit to data allowed a free phase component.

Correlation of inter-glacials with maximum eccentricity of Earth's orbit

Correlation of inter-glacials with maximum eccentricity of Earth

The overall fit to 5.3 million years of proxy temperature data is not perfect but does reproduce the main features of the observed data. The rapid end to recent  Ice Ages resembles better a saw-tooth oscillation rather than a simple harmonic. However the recurrence of Ice Ages over the last 600,000 years is very  well represented by the fit. Therefore this fit can be used with some confidence to extrapolate into the future and thereby predict the onset of the next Ice Age. This is shown below.

Next Ice Age due to start in 5000 years time.

Next Ice Age due to start in 5000 years time.

Detail of near term climate

Detail of near term climate

Within 10,000 years the Earth will be well on the way into a new deep glaciation which should peak 30,000 years into the future. By this time much of North America and Northern Europe will be under a huge ice sheet. Looking in detail at the data it can be seen that already within 2 thousand years from now, the climate will be noticeably cooler than the current climate. It appears to be the case that the climate is at its warmest period and will naturally begin to cool in the future.

Summary:

The earth’s climate has cooled over the last 5 million years leading to a series of major glaciations over the last 900,000 years. Various origins of this cooling have  been proposed such as  the movement of Antarctica over the South Pole and the uplifting of the Tibetan plateau. Detailed analysis of recent  oscillations in temperature show a clear 100,000 year correlation with interglacials coincident with maxima of the elipticity of the Earth’s  orbit. The oscillation is not present more than 1 million years ago, however a 400,000 year super-oscillation of elipticity is evident from 3 – 5.3million years ago. The  41,000 year Milakowitz oscillation is observed throughout the full 5.3 million years with warmer interglacials coincident with a maximum tilt of the the Earth’s axis. No evidence of a significant 23,000 year  oscillation due to  the precession of the axis of rotation is observed.

The parameterized fit to the LR04 data for the last 600,000 years reproduces all previous glaciations and therefore can be used to predict the next Ice Age. It is found that the Earth is currently at its maximum interglacial global temperature and will start cooling within the next 2-5000  years. 20,000 years from now the Earth will likely be in the depths of another major Ice Age.

References

1) Maureen Raymo & Peter Huybers, Unlocking the mysteries of the ice ages, Nature Vol 451/17 P. 284, 2008

2) Concise Dictionary of Physics: Phenomenological Theory. A theory that expresses mathematically the results of observed phenomena without paying detailed attention to their fundamental significance.

3) Lisiecki, L. E., and M. E. Raymo (2005), A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records, Paleoceanography, 20, PA1003

4) Antarctic drift animation http://www.exploratorium.edu/

5) Raymo, M.E., W.F. Ruddiman, and P.N. Froelich (1988) Influence of late Cenozoic mountain building on ocean geochemical cycles. Geology, v. 16, p. 649-653.

6) Gerald H. Haug2 & Ralf Tiedemann, Effect of the formation of the Isthmus of Panama on Atlantic Ocean thermohaline circulation,  Nature 393, 673-676 (18 June 1998)

Posted in Climate Change, Physics | Tagged , | 2 Comments