Climate of the Past Discussions Climate of the Past Discussions is the access reviewed discussion forum of Climate of the Past Astronomical forcing and mathematical theory of glacial-interglacial cycles

There are three important features of a proxy time series recorded during the Late Pleistocene. They are: 1) 100 000-year cycle as a dominant control of global glacialinterglacials through the late Quaternary, 2) fluctuations with periods of about 40 and 20 thousand years (their contribution to dispersion is no more than 20%), 3) “Red5 noise” behavior of the time series. Direct influence of the insolation change created by fluctuations of the eccentricity is too weak to cause the observed 100 000-year climate fluctuations. Therefore, other mechanisms of such a rhythm are proposed. On the basis of the equation of the heat budget, the equation describing dynamics of zonally averaged temperature is developed. Various combinations of terms of this equation 10 are discussed. They present a linear response to the Milankovitch periodicity, the Langeven stochastic equation, the equation of delay oscillator, the stochastic equation of spontaneous transitions, and the equation of stochastic resonance. Orbitally-induced changes in the solar energy flux received by the Earth play an important role as a mechanism starting process of climate changes which is supported 15 and intensified by different feedbacks within the climate system. Positive anomalies of solar radiation serve as a mechanism causing reorganization of the climate only in rare cases when inclination of Earth axis of rotation increases and, simultaneously, perihelion takes place during the summer time (for the northern hemisphere).


Introduction
The climate is determined by the interaction of solar radiation with the atmosphere and terrestrial geography.Long-term orbital variations produce quasi-periodic variations in the seasonal and latitudinal distribution of solar radiation (insolation), while atmospheric composition influences terrestrial heating by controlling the absorption and transmission of incoming insolation and outgoing terrestrial radiation.The configuration of the continents did not change during the Quaternary, apart from the waxing and waning of continental ice sheets and associated changes in sea level and coastal zone geography.Boundary conditions (insolation, greenhouse gas concentrations, land surface changes, and aerosol concentrations) combined to affect atmospheric and oceanic circulation, as well as the extent and height of the ice sheets.Global ice volume lags behind insolation changes by some thousands of years.
One candidate for a forcing agent that could produce cold (glacial) and warm (interglacial) phases is the Milankovitch mechanism (Berger, 1978;Imbrie et al., 1984;Crucifix, 2008).According to this theory, the Earth's orbital parameters change due to the influence of the Moon, Sun, and planets.Eccentricity slowly varies over 100-and 400-kyr periods, inducing small changes in the annual mean total insolation received by the Earth.In spite of this fact, the 100-kyr glacial-interglacial cycle emerging in the geological data is more or less in phase with the eccentricity cycle (high eccentricity at low ice volume, Hays et al., 1976).This reinforces the necessity to call upon other external forcings or leaves open possibility for an internal origin due possibly to nonlinear processes in order to generate this 100-kyr climatic cycle from the astronomical forcing.Obliquity oscillates from 22 • to 25 • over a 41-kyr period, and the position of the equinoxes precesses relative to the perihelion with 19-and 23-kyr periodicities.Obliquity and precession do not lead to global changes in annual mean energy but strongly modulate the seasonal pattern of insolation.Typical spectral density of different core records (both ice and deep-sea) demonstrates that much of the energy at low frequencies corresponds to periods around 100 kyr.At high frequencies (ν), the spectral density shape displays a red-noise curve (∼ν −2 ).The red-noise continuum is superimposed some weak spectral peaks (whose statistical significance needs to be assessed) belonging to the Milankovitch periodicities.On the average, the obliquity band (near 41-kyr period) accounts for less than 11%, and the precessional band (near 20-kyr period) accounts for less than 1% of total variance (Wunsch, 2003).These values of the Milankovitch-band contribution are consistent with, but smaller than, the bound of 25% suggested by Kominz and Pisias (1979).

329
Recent investigations have shown that the 100-kyr periodicity disappears before 1000 kyr ago (Berger et al., 1998;Bintanja and van de Wal, 2008), at a time when ice sheets were much less developed over the Earth.We do not discuss the Mid-Pleistocene 41 to 100-kyr transition and climate conditions during past several million years.Our aim is to discuss the principal link between astronomical parameters and climate changes over the Late Pleistocene when 100-kyr cycle played a dominant role.We plan to introduce simple conceptual mathematical models of climate, developing foundation of theory of climate change during the Quaternary.

Simple mathematical model of climate
We consider a simple model of climate in which variations of the averaged (vertically and zonally) annual-mean temperature T are determined by the radiation balance R at the top of the atmosphere and the atmospheric and oceanic meridional heat transfer.The heat balance equation is given by: where t is the time, ϕ is the latitude, r is the radius of the Earth, F is the meridional heat flux, c 1 and c 2 are constants.The radiation balance is given by: R=0.25I 0 (1−α) −I ↑ , where I 0 is the solar constant, α is the albedo of the atmosphere and surface system, I ↑ is the flux of outgoing (at the top of the atmosphere) longwave terrestrial radiation.Empirical information yields the linear approximation: I ↑ =A+BT .It is assumed that parameterization (1) includes the influence of meridional fluxes of latent heat and potential energy.Such approaches are broadly used in simplified climate models beginning from the energy-balance Budyko-Sellers model.
Consider a temperature variation over a zonally-averaged zone of the hemisphere located between equator and mid-latitude ϕ.Now temperature change is described by ordinary differential equation.It can be solved only when additional expressions are provided for heat fluxes at the boundaries of this area.We assume that the meridional heat flux is given by: F ϕ = − K δT , where δT =T −T ϕ=0 .It is well known that the magnitude of temperature anomalies in the mid-latitude zone is much larger than near the equator.Therefore, we will ignore the contribution of variability of T ϕ=0 to δT .Additionally, let's take into account that the cross-equatorial heat transfer F ϕ=0 =0.As an approximation, the coefficient K is assumed to be given by K ∝ (δT ) 2 , reflecting the fact that the intensity of heat exchange depends on inhomogeneity of the temperature field.Let us subside all expressions into Eq.( 1).Additionally, we can assume that some processes act with a time delay, and corresponding linear (simplest approach) time-delaying term is included into our analyses.Hence, Eq. ( 1) becomes: The last term in Eq. ( 2) symbolically describes the climate response to the perturbation of solar energy received by the Earth due to the Milankovitch effect.The term f is added to this equation in order to take into account a number of uncertainties which accompany the construction of Eqs. ( 1) and ( 2).It is assumed that f is presented by a time series of short-term fluctuations, and their typical period is much less than the timescale of climate change (λ −1 ).Indeed, f is proposed to be a stochastic function described by the model of white noise.
The time-delaying term is used to directly simulate coupled climate variability involving interactions between the components of the Earth system characterized by very different inertia.Indeed, continental ice sheets respond very slowly to temperature signals due to large inertia of huge ice mass, therefore, small influences are integrated, and changes in ice sheets have to be delayed.This behavior can be simulated in different ways, for example, directly by using the delay-coupled oscillator systems.The existence of such and other (global climate -global carbonate system interaction) timedelay processes is supported through observed data.

331
It is clear that the time-delaying expression could be described more correctly in the following form: where σ(0) has a maximum, and this function decreases with increasing s.However, if this decrease is very fast (like the Dirac delta function), instead of Eq. (3), we obtain the expression with fixed time-delay parameter, same as in Eq. ( 2).Analysis of Eq. ( 2) allows us to understand the possible connection between the reconstructed variations and different theoretical models bringing more criteria for accepting or rejecting any proposed link between cause and effect.Relevant information is assumed to be extracted from particular solutions of equations representing various segments of basic Eq. ( 2).Of course, the sum of particular solutions cannot provide general solution of non-linear Eq. ( 2).

Different models of climate change
Let us analyze Eq. (2).First of all, we stress that the linear response of the climate system to the obliquity and precession periodicities involves fluctuations with periods of 41 kyr and both 23 and 19 kyr, respectively.
Let us consider the following combination of the terms of Eq. ( 2): where λ<0 to provide negative feedback.
Equation (4) represents a system with infinite memory and its solution (so-called Wiener process) would be the continuous counterpart of the Brownian motion.Equation ( 5) is the Langeven equation; it depicts the so-called diffusion stochastic process.Such equations were used by Hasselmann (1976) in his stochastic theory of climate and by many later investigators.It is well known that the solution of Eq. ( 5) is a stochastic process, and its spectrum is given by: Hence, at high frequencies (ν λ), the spectral shape displays a red-noise continuum same as depicted by observed data.It appears that much of the climate system dynamics is best understood in terms of such random walks.
The time scale (λ −1 ) has to be several kyr to be comparable with astronomical forcing periodicities.Using an order-of-magnitude analysis of observed paleodata, the time scale can be estimated as λ −1 =T 0 ∆t/∆T .The transition from the last glacial (the late Pleistocene) cold event to the interglacial (the Holocene) lasted ∼5-7 kyr and annual global temperature increased by ∼5 • C (Joussaume, 1999;Kislov et al. 2002).Using T 0 =15 0 C as mean temperature, we obtain λ −1 ≈20 kyr.
Let us consider the following segment of the Eq. ( 2): The construction Ṫ ∝βT 3 , where β<0, plays an important role restricting a rise of anomalies.First, consider simple situation without stochastic forcing (f =0), and Eq. ( 7) is reduced to: It has three equilibrium solutions: T 1,2 =± −λ/β , λ>0, and T 3 =0.Linear analysis shows that T 1 and T 2 are stable states, but T 3 is an unstable state.So, depending on the initial condition, T →T 1 or T →T 2 .

333
If stochastic forcing is strong enough, the solution of Eq. ( 7) has principally another behavior.Even if the solution has practically reached T 1 or T 2 , the stochastic fluctuation can throw T to the neighborhood of another equilibrium state.After that T "is attracted" to a new equilibrium point, but strong forcing can transfer T to the previous state again, and processes are repeated.Therefore, Eq. ( 7) is a model of a so-called "spontaneous transitions".The periods of transitions are random time intervals, but the mean period or frequency transitions (ν tr ) can be evaluated; it depends on the intensity of stochastic forcing.
For numerical calculation, stochastic forcing was approximated by , where ( ) denotes ensemble average value.Its parameters (intensity (c) and scale of correlation (η)) as well λ as β have been chosen so that the solution of the equation has provided the average time of transition ∼50 kyr.It corresponds to the required 100-kyr periodicity.Figure 1 shows the spectral estimates.
Note the contrasting impressions: spectral peak is located near 100 kyr, and the logarithmic scale display gives the impression of a general red-noise process at high frequencies, on which some weak spectral structures are superimposed.An apparent change in the spectrum character occurs across the 100 kyr band, with the process becoming much closer to white noise at longer periods.It is amazing that there is an apparent real enhancement of energy in the vicinity of the obliquity band near one cycle (41 kyr) in spite of the fact that this rhythm is not included into the model (Eq.7).This fact combined with above mentioned information about low contribution of the Milankovitch cycles to full variance, allows us to emphasize that indeed Milankovitch cyclicities are not a prevailing feature of spectra and look as stochastically scattered peaks on a frequency range.
Let us consider the following segment of Eq. ( 2): Solution T is slightly influenced by the stochastic forcing and infinitely increases with increasing t .This term does not produce fluctuations and cannot be used for the simulation of climate cyclisity.Combination of the following terms of Eq. ( 2): (where β<0, λ+α>0) describes the delay oscillator.Its solution (same analysis as for Eq.8) has three equilibrium states: T 1,2 =± −(λ+α)/β and T 3 =0.The problem of stability is investigated using small deviation (T ) from the equilibrium states.Corresponding equation is given by: We look for a solution in the form: T ∝ exp (ξt), determining real and imaginary part of ξ (ξ=ξ r +i ξ i ).The field of values (α, τ) corresponding to unstable regime is bounded by a neutral curve ξ r =0.For any pair of parameters within the unstable zone, solution moves from one equilibrium state to another and visa versa.To numerically simulate the 100-kyr cycle, the delay (τ) is taken as 50 kyr.Steady self-oscillation of temperature is displayed in Fig. 2. All values are normalized using the scale of temperature λ/β.
Let us consider the last segment of Eq. ( 2): The solution of this equation describes the effect of the stochastic resonance providing an example of a noise-induced transition in a nonlinear system driven by an information signal and noise simultaneously.Such nontrivial phenomenon of stochastic synchronization of a bistable system due to external periodic force is realized under conditions when ν k ≈ν tr (Anishchenko et al., 1999).It means that the presence of "spontaneous transitions" within the climate system with ν tr near one cycle/100 kyr (see above) can 335 in fact enhance the weak eccentricity-induced signal of insolation with the same periodicity (Benzi et al., 1982;Nicolis, 1982;Wiesenfeld and Moss, 1995).Finally, it was shown that principal features of observed climate changes can be reproduced by different models.The important role belongs to interaction of elements of climate system with strongly variable inertia.It causes the phenomena look like time-delay processes and the Brownian motion generating such effects as red-noise variability, spontaneous transitions between unstable equilibriums and "cultivation" of weak external signals due to a stochastical resonance.

Note about the Milankovitch effect
The Milankovitch's periodicities provide (as it has been already noted) only less than ∼20% of climate variance at scales of thousand -dozens of thousand years.Formally it means that orbitally-induced insolation changes are primarily important as the mechanism waking up and accelerating internal feedbacks.Therefore this mechanism was defined as "pacemaker" (Hayes et al., 1976) of the global climate changes.Mathematical exercise presented above has demonstrated the mathematical models of this effect.
However, not very often, the climatic response to variations in insolation can be distinguished from noise.It occurs when large obliquity corresponds to a period when the date of perihelion takes place during the Northern hemisphere summer.Received solar energy is increased during summer, and this factor determines a global warming in spite of a negative solar energy anomaly during winter.The opposite situation leads to a global cooling.Hence, in such cases, large insolation anomalies are not only a pacemaker, but they are also an important factor creating essential reorganization of the climate.Climate simulations had demonstrated that this factor was responsible for transition from the cold Late Pleistocene cryochrone to the Holocene (Joussaume, 1999).Another example is transition from the warm marine isotope stage 5e to the cold stage 5d (125-115 kyr BP) (Schlesinger and Verbitsky, 1996).

Conclusions
Much of climate change is describable as a stochastic variable, often associated with random walks.The 100-kyr period is a dominant feature.At periods shorter than 100 kyr, the spectrum becomes red.Weak peaks superimposed upon some of these spectra correspond to the major Milankovitch forcing frequency bands.Reduction of climate system behavior to simple models is obviously a gross oversimplification.Despite this, combination of external forcing and stochastic forcing of a system together with internal non-linear mechanisms (time-delay process, spontaneous transitions between unstable equilibriums and stochastical resonance) can display a variability similar to that observed for the 100 000-year glacial-interglacial oscillation.It is difficult to establish a priority of concrete mechanisms, most likely their joint influence should be considered.
Orbitally-induced insolation changes play a significant role as a pacemaker activating the internal feedbacks of the global climate system and, sometimes, as the main factor determining climate changes.

Fig. 1 .Fig. 2 .
Fig. 1.Power spectrum, S (ν), for the time-series of numerical solution of Eq. (7).A leastsquares fit to the power spectrum over the range of frequencies is indicated by the short dashed line, where S∝ν −2 .At low frequencies, a transition to near white noise occurs.