Orbital variations, insolation, and the ice ages
Contents
Orbital variations, insolation, and the ice ages#
This notebook is part of The Climate Laboratory by Brian E. J. Rose, University at Albany.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from climlab import constants as const
/Users/tompkins/Library/Python/3.10/lib/python/site-packages/climlab/radiation/cam3.py:46: UserWarning: Cannot import and initialize compiled Fortran extension, CAM3 module will not be functional.
warnings.warn('Cannot import and initialize compiled Fortran extension, CAM3 module will not be functional.')
/Users/tompkins/Library/Python/3.10/lib/python/site-packages/climlab/radiation/rrtm/rrtmg_lw.py:19: UserWarning: Cannot import and initialize compiled Fortran extension, RRTMG_LW module will not be functional.
warnings.warn('Cannot import and initialize compiled Fortran extension, RRTMG_LW module will not be functional.')
/Users/tompkins/Library/Python/3.10/lib/python/site-packages/climlab/radiation/rrtm/rrtmg_sw.py:19: UserWarning: Cannot import and initialize compiled Fortran extension, RRTMG_SW module will not be functional.
warnings.warn('Cannot import and initialize compiled Fortran extension, RRTMG_SW module will not be functional.')
/Users/tompkins/Library/Python/3.10/lib/python/site-packages/climlab/convection/emanuel_convection.py:14: UserWarning: Cannot import EmanuelConvection fortran extension, this module will not be functional.
warnings.warn('Cannot import EmanuelConvection fortran extension, this module will not be functional.')
1. The ice ages#
Recent Earth history (past few million years) has been dominated by the repeated growth and retreat of large continental ice sheets, mostly over the land masses of the Northern Hemisphere.
Extent of glaciation#
The images below show typical maximum extents of the ice sheets during recent glaciations (grey) compared with present-day ice sheets (black)
Hannes Grobe/AWI, http://commons.wikimedia.org/wiki/File:Iceage_north-intergl_glac_hg.png
Pacing of ice ages: evidence from ocean sediments#
The figure below shows a global record of oxygen isotopes recorded in the shells of marine organisms. This record tells us primarily about variations in global ice volume – because the net evaporation of water from the oceans to build up the ice sheets leaves the oceans enriched in heavier isotopes.
The x axis is plotted in Thousands of years before present (present-day is at zero on the left).
Lisiecki, L. E. and Raymo, M. E. (2005). A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records. Paleoceanog., 20.
The ice ages (times of extensive glaciation and high ocean \(\delta^{18}\)O) do not seem to be random fluctations. They have come and gone (approximately) periodically, somewhat like the seasons.
Spectral analysis of such records reveals peaks at some special frequencies:
Imbrie, J. and Imbrie, K. P. (1986). Ice Ages: Solving the Mystery. Harvard University Press, Cambridge, Massachusetts.
The peaks noted on this figure are special because they correspond to frequencies of variations in Earth’s orbital parameters, as we will see.
These kind of results became available in the 1970’s for the first time, because ocean sediment cores allowed a sufficiently detailed look into the past to use time series analysis methods on them, e.g. to compute spectra.
The presence of peaks in the spectrum at orbital frequencies was seen as convincing evidence that the so-called astronomical theory of the ice ages was (at least partially) correct.
2. Introducing the astronomical theory of the ice ages#
The Astronomical Theory of climate and the ice ages looks to the regular, predictable variations in the Earth’s orbit around the Sun as the driving force for the growth and melt of the great ice sheets. Such theories have been discussed since long before there was any evidence about the timing of past glaciations.
Last time we saw that insolation is NOT perfectly symmetrically distributed between the two hemispheres and seasons.
To refresh our memory, let’s use
climlab.solar.insolation.daily_insolation()
to compare the maximum insolation received at the North Pole (at its summer solstice) and the South Pole (at its summer solstice).
from climlab.solar.insolation import daily_insolation
days = np.linspace(0, const.days_per_year, 365)
Qnorth = daily_insolation(90,days)
Qsouth = daily_insolation(-90,days)
Qtrieste = daily_insolation(45,days)
Qeq = daily_insolation(0,days)
print( 'Daily average insolation at summer solstice:')
print( 'North Pole: %0.2f W/m2.' %np.max(Qnorth))
print( 'South Pole: %0.2f W/m2.' %np.max(Qsouth))
plt.plot(days,Qnorth,label="North Pole")
plt.plot(days,Qtrieste,label="Trieste")
plt.plot(days,Qeq,label="Equator")
plt.legend()
plt.show()
Daily average insolation at summer solstice:
North Pole: 525.31 W/m2.
South Pole: 562.03 W/m2.
These asymmetries arise because of the detailed shape of the orbit of the Earth around the Sun and the tilt of the Earth’s axis of rotation.
As these orbitals details change over time, there are significant changes in the distribution of sunlight over the seasons and latitudes.
The Milankovitch hypothesis#
Version of the astronomical theory have been debated for at least 150 years.
The most popular flavor has been the so-called Milankovitch hypothesis:
Ice sheets grow during periods of weak summer insolation in the Northern high latitudes.
The idea is that for an ice sheet to grow, seasonal snow must survive through the summer. Milankovitch therefore focussed on the factors determining the climatic conditions during summer.
3. Ellipses and orbits#
First, watch this neat animation from Peter Huybers (Harvard University):
http://www.people.fas.harvard.edu/~phuybers/Inso/Orbit.mv4
Watch carefully and note the three ways that the orbit is varying simultaneously.
From Professor Huybers’ web page:
A movie depicting Earth’s changing orbit over the last 100Ky. The orientation is such that spring equinox (indicated by a vertical bar) is directly to the front with the sun behind it. Northern Hemisphere summer is to our right, and Northern Hemisphere winter is to the left. The apsidal (dashed) line connects perihelion (Earth’s closest approach to the sun) to aphelion (the point when Earth is furthest from the sun). The rotaion of the apsidal line occurs because of the precession of the equinoxes and has a roughly twenty-two thousand year period. The semi-circle around the Earth indicates the location of the equator and the straight line is the polar axis. Obliquity is defined as the angle beetween the orbital and equatorial planes. The variations in Earth’s obliquity and the eccentricity of Earth’s orbit have both been increased in magnitude by a factor of ten. Also, the Earth’s angular velocity has been decreased by a factor of five thousand. Note that Earth’s angular velocity is slowest at aphelion and fastest at perihelion.
The Earth’s orbit around the Sun traces out an ellipse, with the Sun at one focal point.
Perihelion and Aphelion#
The point in the orbit that is closest to the sun is called Perihelion. The farthest point is called Aphelion.
Distances (present-day):
Perihelion, \( d_p = 1.47 \times 10^{11}\) m
Aphelion, \( d_a = 1.52 \times 10^{11}\) m
Eccentricity#
The eccentricity of the orbit is defined as $\( e = \frac{d_a-d_p}{d_a+d_p} \)$
So for present-day values, \(e = 0.017 = 1.7\%\)
Earth’s orbit is nearly circular, but not quite!
(What value of \(e\) would a purely circular orbit have?)
As the Earth travels around its orbit, the distance to the sun varies. The energy flux (W m\(^{-2}\)) is larger when the Earth is closer to the sun (i.e. near perihelion).
At present, perihelion occurs on January 3. This is very close to the Northern Hemisphere winter solstice (Dec. 21).
The Earth actually receives MORE total sunlight during Northern Winter than during Northern Summer.
It is thus critical to understand the relative timing of our seasons (which are determined by the axial tilt or obliquity) and the perehelion.
Obliquity#
The obliquity \(\Phi\) is the tilt of the Earth’s axis of rotation with respect to a line normal to the plane of the Earth’s orbit around the Sun.
Currently \(\Phi = 23.5^\circ\)
Obliquity is the fundamental reason we have seasons, and would have seasons even with a perfectly circular orbit (\(e=0\))
Higher obliquity means:
more summertime insolation at the poles
less wintertime insolation in mid-latitudes
Longitude of perihelion and precession of the equinoxes#
The longitude of perihelion is defined as the angle \(\Lambda\) between the Earth-Sun line at vernal equinox and the line from the Sun to perihelion (see sketch).
The current value is \(\Lambda = 281^\circ\) (perihelion on January 3, shortly after NH winter solstice).
We call the gradual change over time of the longitude of perihelion the precession of the equinoxes (or just precession). It is the gradual change in the time of year at which the Earth is closest to the Sun.
It is important to understand that eccentricity modulates the precession. Highly eccentric orbits lead to larger differences in the seasonal distribution of insolation.
We quantify this with the precessional parameter
Large positive precessional parameter = Excess insolation during summer in the northern hemisphere.
The three orbital parameters#
We have just identified three parameters that control the seasonal and latitudinal distribution of insolation: \(e, \Lambda, \Phi\)
All three vary in predictable ways over time. They have been calculated very accurately from astronomical considerations (basically the gravity of the Earth, Sun, moon, and other solar system objects).
4. Past orbital variations#
There are tools in climlab
to look up orbital parameters for Earth over the last 5 million years.
We will use the package
climlab.solar.orbital
from climlab.solar.orbital import OrbitalTable
Tokenization took: 3.31 ms
Type conversion took: 2.82 ms
Parser memory cleanup took: 0.00 ms
OrbitalTable
<xarray.Dataset> Dimensions: (kyear: 5001) Coordinates: * kyear (kyear) int64 0 -1 -2 -3 -4 -5 ... -4996 -4997 -4998 -4999 -5000 Data variables: ecc (kyear) float64 0.01724 0.01764 0.01802 ... 0.02689 0.02685 long_peri (kyear) float64 281.4 264.3 247.2 ... -8.361e+04 -8.362e+04 obliquity (kyear) float64 23.45 23.57 23.7 23.82 ... 22.55 22.52 22.52 precession (kyear) float64 -0.0169 -0.01756 -0.01662 ... -0.02683 -0.02639 65NJul (kyear) float64 426.8 430.1 434.7 440.2 ... 407.5 409.1 412.3 65SJan (kyear) float64 455.6 455.5 454.1 451.6 ... 452.9 450.6 447.3 15NJul (kyear) float64 440.6 442.5 445.6 449.8 ... 431.5 433.5 436.9 15SJan (kyear) float64 470.4 468.6 465.5 461.4 ... 479.7 477.5 474.0 Attributes: Description: The Berger and Loutre (1991) orbital data table Citation: https://doi.org/10.1016/0277-3791(91)90033-Q Source: http://www.atmos.albany.edu/facstaff/brose/resources/climla... Note: Longitude of perihelion is defined to be 0 degrees at North...
Make reference plots of the variation in the three orbital parameter over the last 1 million years
kyears = np.arange( -1000., 1.)
print(kyears)
orb = OrbitalTable.interp(kyear=kyears)
orb
[-1000. -999. -998. ... -2. -1. 0.]
<xarray.Dataset> Dimensions: (kyear: 1001) Coordinates: * kyear (kyear) float64 -1e+03 -999.0 -998.0 -997.0 ... -2.0 -1.0 0.0 Data variables: ecc (kyear) float64 0.03576 0.03695 0.03811 ... 0.01764 0.01724 long_peri (kyear) float64 -1.644e+04 -1.642e+04 -1.641e+04 ... 264.3 281.4 obliquity (kyear) float64 23.78 23.84 23.88 23.9 ... 23.7 23.57 23.45 precession (kyear) float64 0.03018 0.02458 0.0166 ... -0.01756 -0.0169 65NJul (kyear) float64 478.7 479.0 476.4 470.9 ... 434.7 430.1 426.8 65SJan (kyear) float64 414.9 416.2 419.9 425.8 ... 454.1 455.5 455.6 15NJul (kyear) float64 489.6 489.1 485.9 480.0 ... 445.6 442.5 440.6 15SJan (kyear) float64 424.4 425.0 428.3 434.0 ... 465.5 468.6 470.4 Attributes: Description: The Berger and Loutre (1991) orbital data table Citation: https://doi.org/10.1016/0277-3791(91)90033-Q Source: http://www.atmos.albany.edu/facstaff/brose/resources/climla... Note: Longitude of perihelion is defined to be 0 degrees at North...
The object called
orb
is an xarray.Dataset
which now holds 1 million years worth of orbital data, total of 1001 data points for each element:
eccentricity
ecc
obliquity angle
obliquity
solar longitude of perihelion
long_peri
fig = plt.figure( figsize = (10,12) )
ax1=fig.add_subplot(3,1,1)
ax1.plot( kyears, orb['ecc'] )
ax1.set_title('Eccentricity $e$', fontsize=18)
ax2 = fig.add_subplot(3,1,2)
ax2.plot( kyears, orb['ecc'] * np.sin( np.deg2rad( orb['long_peri'] ) ) )
ax2.set_title('Precessional parameter $e \sin(\Lambda)$', fontsize=18 )
ax3 = fig.add_subplot(3,1,3)
ax3.plot( kyears, orb['obliquity'] )
ax3.set_title('Obliquity (axial tilt) $\Phi$', fontsize=18 )
ax3.set_xlabel( 'Thousands of years before present', fontsize=14 )
plt.savefig('plot.png')
Timescales of orbital variation:#
Eccentricity varies slowly between nearly circular and slightly eccentric, with dominant periodicities of about 100 and 400 kyear. Current eccentricity is relatively small compared to previous few million years
Longitude of perihelion has a period around 20 kyears, but effect is modulated by slow eccentricity variations. Precessional cycles are predicted to be small for the coming 50 kyears because of weak eccentricity.
Obliquity varies between about 22.5º and 24.5º over a period of 40 kyears. It is currently near the middle of its range.
5. Using climlab
to calculate insolation for arbitrary orbital parameters#
We can use the function
climlab.solar.insolation.daily_insolation()
to calculate insolation for any arbitrary orbital parameters.
We just need to pass a dictionary of orbital parameters. This works automatically when we slice or interpolate from the xarray
object OrbitalTable
.
An example: zero obliquity#
Calculate the insolation at the North Pole for a planet with zero obliquity and eccentricity.
from climlab.solar.insolation import daily_insolation
thisorb = {'ecc':0.0, 'obliquity':0., 'long_peri':0.}
print(thisorb)
{'ecc': 0.0, 'obliquity': 0.0, 'long_peri': 0.0}
days = np.linspace(1.,20.)/20 * const.days_per_year
daily_insolation(90, days, thisorb)
<xarray.DataArray (day: 50)> array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) Coordinates: * day (day) float64 18.26 25.34 32.42 39.51 ... 344.0 351.1 358.2 365.2
Compare this with the same calculation for default (present-day) orbital parameters:
daily_insolation(90, days)
<xarray.DataArray (day: 50)> array([ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 18.86214697, 85.00612352, 149.06416099, 210.14318852, 267.43771123, 320.2345592 , 367.91459762, 409.95198005, 445.91155305, 475.44500661, 498.28632029, 514.24699096, 523.21145124, 525.13300552, 520.03052731, 507.9860795 , 489.14354007, 463.70823983, 431.94754511, 394.19224658, 350.83854561, 302.35036042, 249.26160895, 192.17806354, 131.7783206 , 68.81338909, 4.10438408, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ]) Coordinates: * day (day) float64 18.26 25.34 32.42 39.51 ... 344.0 351.1 358.2 365.2
Do you understand what’s going on here?
6. Past changes in insolation: investigating the Milankovitch hypothesis#
The Last Glacial Maximum or “LGM” occurred around 23,000 years before present, when the ice sheets were at their greatest extent.
By 10,000 years ago, the ice sheets were mostly gone and the last ice age was over.
If the Milankovitch hypothesis is correct, we should that summer insolation in the high northern latitudes increased substantially after the LGM.
The classical way to plot this is the look at insolation at summer solstice at 65ºN. Let’s plot this for the last 100,000 years.
# Plot summer solstice insolation at 65ºN
kyears = np.linspace(-100, 0, 101) # last 100 kyr
thisorb = OrbitalTable.interp(kyear=kyears)
# first 65=latitude
# second 172 day of year
S65 = daily_insolation( 65, 172, thisorb )
print(S65)
fig, ax = plt.subplots()
ax.plot(kyears, S65)
ax.set_xlabel('Thousands of years before present')
ax.set_ylabel('W/m2')
ax.set_title('Summer solstice insolation at 65N')
ax.grid()
<xarray.DataArray (kyear: 101)>
array([498.80375813, 489.01652477, 480.9736246 , 475.3502372 ,
472.2993669 , 471.77014222, 473.41491328, 476.88238506,
481.71637865, 487.5535188 , 494.02933045, 500.8172215 ,
507.60502296, 514.03794513, 519.72839453, 524.18459749,
526.88818882, 527.31498761, 525.12731424, 520.26440099,
512.97325288, 503.90310504, 493.88961057, 483.85946184,
474.60972819, 466.80827497, 460.84556162, 456.91158754,
455.01008584, 455.05498997, 456.83533549, 460.10625584,
464.61664832, 470.0812763 , 476.24469622, 482.76915403,
489.37529882, 495.71509717, 501.47781924, 506.35086965,
510.11431306, 512.6118726 , 513.84752665, 513.88974921,
512.96185562, 511.29158609, 509.19938758, 506.92430447,
504.68831494, 502.63129842, 500.87843695, 499.44012914,
498.31547751, 497.5029162 , 496.96147228, 496.63872361,
496.49184725, 496.50663619, 496.64846569, 496.87392298,
497.17828287, 497.49832129, 497.76826152, 497.90450296,
497.79616493, 497.28260953, 496.21327019, 494.43784106,
491.84687581, 488.44731105, 484.33284194, 479.78139543,
475.13202835, 470.82223114, 467.30625982, 464.94659708,
464.04012809, 464.71360749, 467.01182992, 470.83601977,
476.0017198 , 482.28701237, 489.37366141, 496.92573779,
504.59665857, 511.92465686, 518.48117328, 523.7717273 ,
527.3444343 , 528.82231027, 527.97485582, 524.87702832,
519.81679294, 513.38066041, 506.22510927, 499.05977296,
492.51031214, 487.03935919, 482.90505474, 480.2196011 ,
478.94375843])
Coordinates:
* kyear (kyear) float64 -100.0 -99.0 -98.0 -97.0 ... -3.0 -2.0 -1.0 0.0
Indeed, there was an increase of 60 W m\(^{-2}\) over a 10 kyr interval following the LGM.
Why?
What orbital factors favor high insolation at 65ºN at summer solstice?
high obliquity (greater tilt towards sun)
large, positive precessional parameter (Earth closest to sun)
Looking back at our plots of the orbital parameters, it turns out that both were optimal around 10,000 years ago.
Actually 10,000 years ago the climate was slightly warmer than today and the ice sheets had mostly disappeared already.
The LGM occurred near a minimum in summer insolation in the north – mostly due to obliquity reaching a minimum, since we have been in a period of weak precession due to the nearly circular orbit. So this is consistent with the orbital theory.
The hypothesis is incomplete, but compelling.
Comparing insolation at 10 kyr and 23 kyr#
lat = np.linspace(-90, 90, 181)
days = np.linspace(1.,50.)/50 * const.days_per_year
orb_0 = OrbitalTable.interp(kyear=0) # present-day orbital parameters
orb_10 = OrbitalTable.interp(kyear=-10) # orbital parameters for 10 kyrs before present
orb_23 = OrbitalTable.interp(kyear=-23) # 23 kyrs before present
Q_0 = daily_insolation( lat, days, orb_0 )
Q_10 = daily_insolation( lat, days, orb_10 ) # insolation arrays for each of the three sets of orbital parameters
Q_23 = daily_insolation( lat, days, orb_23 )
# Qlist=[]
# for y in [0,-10,-23]:
# Qlist.append(daily_insolation(lat,days,OrbitalTable.interp(kyear=y)))
fig = plt.figure( figsize=(12,6) )
ax1 = fig.add_subplot(1,2,1)
Qdiff = Q_10 - Q_23
CS1 = ax1.contour( days, lat, Qdiff, levels = np.arange(-100., 100., 10.) )
ax1.clabel(CS1, CS1.levels, inline=True, fmt='%r', fontsize=10)
ax1.contour( days, lat, Qdiff, levels = [0], colors='k' )
ax1.set_xlabel('Days since January 1', fontsize=16 )
ax1.set_ylabel('Latitude', fontsize=16 )
ax1.set_title('Insolation differences: 10 kyrs - 23 kyrs', fontsize=18 )
ax2 = fig.add_subplot(1,2,2)
ax2.plot( np.mean( Qdiff, axis=1 ), lat )
ax2.set_xlabel('W m$^{-2}$', fontsize=16 )
ax2.set_ylabel( 'Latitude', fontsize=16 )
ax2.set_title(' Annual mean differences', fontsize=18 )
ax2.set_ylim((-90,90))
ax2.grid()
This figure shows that the insolation at summer solstice does not tell the whole story!
For example, the insolation in late summer / early fall apparently got weaker between 23 and 10 kyr (in the high northern latitudes).
The annual mean plot is perfectly symmetric about the equator.
This actually shows a classic obliquity signal: at 10 kyrs, the axis close to its maximum tilt, around 24.2º. At 23 kyrs, the tilt was much weaker, only about 22.7º. In the annual mean, a stronger tilt means more sunlight to the poles and less to the equator. This is very helpful if you are trying to melt an ice sheet.
Finally, take the global average of the difference. Recall that the cells reduce in size as we move towards the poles, so if we take the simple arithmetic mean we will get the incorrect answer, as it will bias the result towards the poles.
print( np.average(np.mean(Qdiff,axis=1), weights=np.cos(np.deg2rad(lat))) )
0.006510430783267399
The difference is tiny (and due to very small changes in the eccentricity).
Ice ages are driven by seasonal and latitudinal redistributions of solar energy, NOT by changes in the total global amount of solar energy!
EXERCISE#
If you recall, 6 Kyr ago, the Sahara was covered by vegetation, and the monsoon penetrated far further north over Africa. Repeat the previous exercise but comparing the isolation 6000 years ago to the present day.
this is a note?
How much more insolation was there at 40N at the summer solstice?
7. Understanding the effects of orbital variations on insolation#
We are going to create a figure showing past time variations in three quantities:
Global, annual mean insolation
Annual mean insolation at high northern latitudes
Summer solstice insolation at high northern latitudes
which we will compare to the orbital variations we plotted earlier.
Create a large array of insolation over the whole globe, whole year, and for every set of orbital parameters.
print(days)
[ 7.304844 14.609688 21.914532 29.219376 36.52422 43.829064
51.133908 58.438752 65.743596 73.04844 80.353284 87.658128
94.962972 102.267816 109.57266 116.877504 124.182348 131.487192
138.792036 146.09688 153.401724 160.706568 168.011412 175.316256
182.6211 189.925944 197.230788 204.535632 211.840476 219.14532
226.450164 233.755008 241.059852 248.364696 255.66954 262.974384
270.279228 277.584072 284.888916 292.19376 299.498604 306.803448
314.108292 321.413136 328.71798 336.022824 343.327668 350.632512
357.937356 365.2422 ]
lat = np.linspace(-90, 90, 91)
num = 365.
days = np.linspace(1.,num,365)/num * const.days_per_year
Q = daily_insolation(lat, days, orb)
print( Q.shape)
(91, 365, 1001)
Qann = np.mean(Q, axis=1) # time average over the year
print(Qann.shape)
print (kyears)
Qglobal = np.empty_like( kyears )
for n in range(kyears.shape[0] ): # global area-weighted average
Qglobal[n] = np.average( Qann[:,n], weights=np.cos(np.deg2rad(lat)))
print( Qglobal.shape)
(91, 1001)
[-100. -99. -98. -97. -96. -95. -94. -93. -92. -91. -90. -89.
-88. -87. -86. -85. -84. -83. -82. -81. -80. -79. -78. -77.
-76. -75. -74. -73. -72. -71. -70. -69. -68. -67. -66. -65.
-64. -63. -62. -61. -60. -59. -58. -57. -56. -55. -54. -53.
-52. -51. -50. -49. -48. -47. -46. -45. -44. -43. -42. -41.
-40. -39. -38. -37. -36. -35. -34. -33. -32. -31. -30. -29.
-28. -27. -26. -25. -24. -23. -22. -21. -20. -19. -18. -17.
-16. -15. -14. -13. -12. -11. -10. -9. -8. -7. -6. -5.
-4. -3. -2. -1. 0.]
(101,)
print(orb)
<xarray.Dataset>
Dimensions: (kyear: 1001)
Coordinates:
* kyear (kyear) float64 -1e+03 -999.0 -998.0 -997.0 ... -2.0 -1.0 0.0
Data variables:
ecc (kyear) float64 0.03576 0.03695 0.03811 ... 0.01764 0.01724
long_peri (kyear) float64 -1.644e+04 -1.642e+04 -1.641e+04 ... 264.3 281.4
obliquity (kyear) float64 23.78 23.84 23.88 23.9 ... 23.7 23.57 23.45
precession (kyear) float64 0.03018 0.02458 0.0166 ... -0.01756 -0.0169
65NJul (kyear) float64 478.7 479.0 476.4 470.9 ... 434.7 430.1 426.8
65SJan (kyear) float64 414.9 416.2 419.9 425.8 ... 454.1 455.5 455.6
15NJul (kyear) float64 489.6 489.1 485.9 480.0 ... 445.6 442.5 440.6
15SJan (kyear) float64 424.4 425.0 428.3 434.0 ... 465.5 468.6 470.4
Attributes:
Description: The Berger and Loutre (1991) orbital data table
Citation: https://doi.org/10.1016/0277-3791(91)90033-Q
Source: http://www.atmos.albany.edu/facstaff/brose/resources/climla...
Note: Longitude of perihelion is defined to be 0 degrees at North...
fig = plt.figure(figsize = (16,12))
ax = []
for n in range(6):
ax.append(fig.add_subplot(3,2,n+1))
ax[0].plot( kyears, orb['ecc'] )
ax[0].set_title('Eccentricity $e$', fontsize=18 )
ax[2].plot( kyears, orb['obliquity'] )
ax[2].set_title('Obliquity (axial tilt) $\Phi$', fontsize=18 )
ax[4].plot( kyears, orb['ecc'] * np.sin( np.deg2rad( orb['long_peri'] ) ) )
ax[4].set_title('Precessional parameter $e \sin(\Lambda)$', fontsize=18 )
ax[1].plot( kyears, Qglobal )
ax[1].set_title('Global, annual mean insolation', fontsize=18 )
ax[1].ticklabel_format( useOffset=False )
ax[3].plot( kyears, Qann[80,:] )
ax[3].set_title('Annual mean insolation at 70N', fontsize=18 )
ax[5].plot( kyears, Q[80,170,:] )
ax[5].set_title('Summer solstice insolation at 70N', fontsize=18 )
for n in range(6):
ax[n].grid()
for n in [4,5]:
ax[n].set_xlabel( 'Thousands of years before present', fontsize=14 )
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[19], line 6
3 for n in range(6):
4 ax.append(fig.add_subplot(3,2,n+1))
----> 6 ax[0].plot( kyears, orb['ecc'] )
7 ax[0].set_title('Eccentricity $e$', fontsize=18 )
8 ax[2].plot( kyears, orb['obliquity'] )
File ~/Library/Python/3.10/lib/python/site-packages/matplotlib/axes/_axes.py:1662, in Axes.plot(self, scalex, scaley, data, *args, **kwargs)
1419 """
1420 Plot y versus x as lines and/or markers.
1421
(...)
1659 (``'green'``) or hex strings (``'#008000'``).
1660 """
1661 kwargs = cbook.normalize_kwargs(kwargs, mlines.Line2D)
-> 1662 lines = [*self._get_lines(*args, data=data, **kwargs)]
1663 for line in lines:
1664 self.add_line(line)
File ~/Library/Python/3.10/lib/python/site-packages/matplotlib/axes/_base.py:311, in _process_plot_var_args.__call__(self, data, *args, **kwargs)
309 this += args[0],
310 args = args[1:]
--> 311 yield from self._plot_args(
312 this, kwargs, ambiguous_fmt_datakey=ambiguous_fmt_datakey)
File ~/Library/Python/3.10/lib/python/site-packages/matplotlib/axes/_base.py:504, in _process_plot_var_args._plot_args(self, tup, kwargs, return_kwargs, ambiguous_fmt_datakey)
501 self.axes.yaxis.update_units(y)
503 if x.shape[0] != y.shape[0]:
--> 504 raise ValueError(f"x and y must have same first dimension, but "
505 f"have shapes {x.shape} and {y.shape}")
506 if x.ndim > 2 or y.ndim > 2:
507 raise ValueError(f"x and y can be no greater than 2D, but have "
508 f"shapes {x.shape} and {y.shape}")
ValueError: x and y must have same first dimension, but have shapes (101,) and (1001,)
We see that
Global annual mean insolation varies only with eccentricity (slow), and the variations are very small!
Annual mean insolation varies with obliquity (medium). Annual mean insolation does NOT depend on precession!
Summer solstice insolation at high northern latitudes is affected by both precession and obliquity. The variations are large.
8. Summary#
The annual, global mean insolation varies only as a result of eccentricity \(e\). The changes are very small (about 0.1 % through a typical eccentricity cycle from more circular to more elliptical)
Obliquity controls the annual-mean equator-to-pole insolation gradient.
The precessional parameter \(e \sin\Lambda\) controls the modulation in seasonal insolation due to eccentricity and longitude of perihelion \(\Lambda\).
The combined effects can result in 15% changes in high-latitude summer insolation
Obliquity combined with eccentricity and longitude of perihelion control the amplitude of seasonal insolation variations at a point.
Combined effects of the three orbital parameters can cause variations in seasonal insolation as large as 30% in high latitudes.
These geometrical considerations tell us that seasonal variations in \(Q\) can be rather large, and will surely impact the climate. But to go from there to understanding how large ice sheets come and go is a difficult step, and requires climate models!
One thing is clear: any serious astronomical theory of climate needs to take account of the climate response to seasonal variations in Q, because these are much larger than the variations in annual mean insolation.
Credits#
Adrian Tompkins forked this notebook April 23rd 2021.
This notebook is part of The Climate Laboratory, an open-source textbook developed and maintained by Brian E. J. Rose, University at Albany.
It is licensed for free and open consumption under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.
Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.