1 Introduction

A major factor in the pacing of climate change on time scales of tens of thousands of years seems to be the change in the seasonal sunlight distribution induced by variations in the Earth’s orbital parameters (e.g., Milankovitch 1941; Berger 1978). The shape of the Earth’s orbit shifts from more elliptical to more circular at periodicities of about 100 ky (ky, 1000 years). Other rhythms are linked to the tilt of the Earth’s axis (periodicities ∼40 ky) and the timing of the seasons relative to the Earth’s closest position to the Sun on the orbit (periodicities ∼20 ky).

Various terrestrial and marine records show that the Holocene and the last interglacial were marked by coherent patterns of climate variability at a regional or global scale (e.g., Keigwin and Pickart 1999; Rimbu et al. 2003a). Pollen-based temperature reconstructions suggest a temperature drop from 120 to 115 kyBP (kyBP, 1000 years before present) in northern and central Europe for boreal winter and summer (Kühl and Litt 2003). Coral data from the northern Red Sea suggest a modulation of the seasonal cycle for the Eemian and late Holocene climate (Felis et al. 2004). These patterns of the past climate variability are often related to changes in major phenomena that dominated the climate variability during the last century, like the El Niño Southern Oscillation (ENSO) (Clement et al. 1999; Tudhope et al. 2001; Kitoh and Murakami 2002), the Arctic Oscillation/North Atlantic Oscillation (AO/NAO) (Glowienka 1985; Hurrell 1995; Thompson and Wallace 1998), the Pacific/North American (PNA) pattern (Wallace and Gutzler 1981), and the Northern/Southern annular modes (Hartmann and Lo 1998; Thompson and Wallace 1998, 2000; Thompson et al. 2000).

We address the question whether these modes of variability are relevant on orbital time scales for the last 140,000 years. The spatial heterogeneity of the interglacial climate evolution is examined by applying a general circulation model and addressing the modulation of large-scale climate modes on orbital time scales.

2 Methodology

For the simulation of interglacial climate changes owing to orbital forcing, we apply the coupled atmosphere–ocean general circulation model ECHO-G (Legutke and Voss 1999). The atmospheric part of this model is the circulation model ECHAM4 with T30 horizontal resolution (approximately 3.8∘ × 3.8∘) and 19 vertical levels. The ocean model includes a dynamic-thermodynamic sea-ice model with snow cover and has a resolution of approximately 2.8∘ × 2.8∘. In the tropics, its meridional resolution is increased to 0.5 ∘. The model consists of 20 irregularly spaced vertical levels with 10 levels covering the upper 300 m.

The control experiment has been integrated over 3000 years of model simulation into a pre-industrial climate state, which is regarded as the quasi-equilibrium response of the model to interglacial boundary conditions (Lorenz and Lohmann 2004). It serves as a baseline and initialization of the simulations where the influence of the solar radiation due to varying orbital parameters has taken into account (Berger 1978). Figure 1 shows the solar insolation variation at different latitudes and times of the year. The marine isotope stages (MIS) are alternating warm and cool periods deduced from oxygen isotope data (Hays et al. 1976). MIS 5 is a warm interglacial (with subperiods of 5a-5a).

Fig. 1
figure 1

Solar insolation variation at different latitudes and times of the year due to varying orbital parameters (Berger 1978). Eq. stands for equator. The insolation variation is the only forcing in the model (see text). The numbers at the top indicate the common numbers to label the marine isotope stages (MIS)

Full size image

The experiment covers the time period of the last 140,000 years to the next 30,000 years into the future. Computer resources for running a complex model like ECHO-G over 170 ky are very demanding. Therefore, the time scale of the orbital forcing has been shortened by an acceleration factor of 100, which implies that the response of the system to the insolation forcing is completed after decades. The 170 ky are therefore represented in 1700 model years. In effect, this procedure decouples any slowly responding portions of the system (e.g., deep water properties) from the feedback structure (Lorenz and Lohmann 2004). Throughout the whole experiments, fixed pre-industrial greenhouse gas concentrations and modern values for vegetation, continental ice distribution, and sea level were used.

Finally, we will compare the simulated temperature evolution with the Eemain and Holocene (MIS 5e and MIS 1 in Fig. 1). Because the dating uncertainties for most marine proxy data are difficult when going further back than 60,000 years (e.g., Drysdale et al. 2009), we will concentrate on the Holocene for the data-model comparison. The proxy data were focused on updating of the global database for proxy-derived Holocene sea surface temperature (SST) records, i.e., an SST synthesis based on alkenone-derived SST estimation and a new synthesis effort compiling the SST records derived from foraminiferal Mg/Ca (Leduc et al. 2010).

3 Results

Since we are interested in the atmospheric dynamics, we concentrate on boreal winter circulation changes (Fig. 2). An AO/NAO index is calculated from the sea level pressure (SLP) region located at Iceland. The results do not sensitively depend on the chosen region of the index. In order to suppress the noise in the system, the first 10 vectors of a Singular Spectrum Analysis (Ghil et al. 2002) of the index were taken (lag; 100 years, the results are independent of the exact value). In our case, this procedure is equivalent to a low pass filter of the signal. The time series is highly negatively correlated with the precipitation in the tropical Pacific Ocean (Fig. 2) which is a sensitive measure of the convection and strength of the tropical atmospheric circulation.

Fig. 2
figure 2

Normalized index time series for boreal winter (DJF). Black line indicates AO/NAO index based on the SLP for the region (30 ∘–0 ∘ W; 55 ∘–80 ∘ N). Red line indicates precipitation index for the region (150 ∘ W–150 ∘ E; 5 ∘ S–5 ∘ N)

Full size image

We apply composite mapping to display a climate sensitivity with respect to these indices. The data were separated as occurring during a high index year, a negative index year, or a neutral year (when it is between ± 1 standard deviation). The data within these three categories were then averaged, and the differences between the positive and negative average were taken for each grid point. The pictures depict an average anomaly field with respect to the index time series. The results for the precipitation index are displayed, the resulting composite maps with respect to the other indices in Fig. 2 are almost identical.

Figure 3 displays the composite map of the boreal winter 500 hPa geopotential height field related to the normalized precipitation index in the tropical Pacific Ocean shown in Fig. 2. The convection anomaly and its associated diabatic heating directly affects the Hadley circulation and the Northern Annular Mode (NAM) pattern. We find that persistent tropical positive (negative) anomalies induce a poleward (equatorward) displacement of the mid-latitude storm tracks (Fig. 3). The latitudinal displacement of the Northern Hemisphere storm track is linked with a negative AO/NAO and NAM phase. The global picture of the 500 hPa geopotential height composite field indicates a global response to low-latitude forcing and bears similarity with the negative phase of Antarctic Oscillation (AAO) in the Southern Hemisphere. The AAO is the dominant pattern of non-seasonal tropospheric circulation variations south of 20 ∘ S (Marshall 2002; Thompson and Wallace 2000). The AAO is also referred to as the Southern Annular Mode (SAM).

Fig. 3
figure 3

Composite map of the boreal winter 500 hPa geopotential height field related to normalized precipitation index in the tropical Pacific Ocean (shown as red curve in Fig. 2). Upper panel indicates Northern Hemisphere polar stereographic projection emphasizing the negative phase of NAO (NAO- pattern). Lower panel indicates global projection indicating global response to tropical forcing. Units are m. Values above 10 and below −18 m exceed the 95 % significance level of a Students’s t test

Full size image

Figure 4 shows the composite maps of SLP field and surface temperature related to the normalized tropical precipitation index. The boreal winter composite map depicts a spatial signature in the SLP and temperature fields. The temperature panel indicates El Niño conditions. The figure reveals a modulation on orbital time scales, mainly driven by precession at low latitudes. When calculating local temperature evolutions based on Fig. 4, we find dominant ∼20 ky variations linked to precession for regions sensitive to the SLP pattern, whereas for the zonal mean (60 ∘–90 ∘ N) the temperature variations show also the strong ∼40 ky obliquity cycle.

Fig. 4
figure 4

Composite map of the boreal winter SLP (upper panel) and surface temperature (lower panel) related to normalized precipitation index in the tropical Pacific Ocean shown in Fig. 2. Units are hPa and ∘C per normalized precipitation index. Values above 1 and below −1 hPa and above 0.5 and below −0.5 ∘C exceed the 95 % significance level of a Students’ t test

Full size image

Figure 5 indicates the composite map of annual mean surface temperature related to the normalized precipitation index in the tropical Pacific Ocean. It clearly shows the NAM-/NAO- signature, seen in boreal winter (Fig. 4). The related atmospheric circulation patterns induce non-uniform temperature anomalies, much stronger in amplitude than by the direct solar insolation.

Fig. 5
figure 5

Composite map of annual mean surface temperature related to the normalized precipitation index in the tropical Pacific Ocean shown in Fig 2. Units are ∘C per normalized precipitation index

Full size image

It is useful to look for the mechanisms of glacial inception at the end of the last interglacial (MIS 5e), as simulated by the ECHO-G model. Figure 6 shows the boreal winter and summer anomalies for 115 minus 120 kyBP. In concert with the response to direct summer solar insolation (Fig. 6c), the atmospheric circulation pattern with the negative phase of NAO (Fig. 6b) provides for a temperature drop over the Northern Eurasian continent and Arctic Ocean for 115,000 years before present (Fig. 6a). This results in an annual mean cooling (Fig. 6d).

Fig. 6
figure 6

Surface air temperature for 115 minus 120 kyBP, for winter (DJF), summer (JJA), and annual mean. For boreal winter (DJF), the sea level pressure anomaly is shown (panel b). All values represent an average over 20 years of the simulation (representing 2000 years) centered at 115 and 120 kyBP, respectively. Units are ∘C

Full size image

A schematic diagram of northern hemisphere forcing (Fig. 7) summarizes our results: A high insolation at low latitudes (blue line in Fig. 1) provides a cooling for northern Eurasia due to atmospheric dynamics during boreal winter (Figs. 3, 4, and 6a, b). Low insolation during boreal summer (green line in Fig. 1) provides a summer cooling (Fig. 6c). Both effects are contributing to an annual mean cooling as seen in Figs. 5 and 6d.

Fig. 7
figure 7

Schematic diagramme of Northern Hemisphere forcing. Left indicates a high insolation at low latitudes (blue line in Fig. 1) provides a cooling for northern Eurasia due to atmospheric dynamics during boreal winter (Figs. 3, 4, and 6a, b). Right indicates a low insolation during boreal summer (green line in Fig. 1) provides a summer cooling (Fig. 6c). Both effects are contributing to an annual mean cooling as seen in Figs. 5, 6d

Full size image

In order to compare the model result with paleoclimate data, we compare furthermore the annual mean SST trends from the mid-Holocene to the present (6000 kyBP) in the simulations and as estimated from alkenone and Mg/Ca temperature proxies for the same time period (Fig. 8). The marine alkenone-based temperature reconstructions (Leduc et al. 2010) comprise marine proxy records for SST based on alkenones and Mg/Ca. The data are unevenly distributed over the world ocean and are mainly located in the North Atlantic Ocean and in coastal areas. Lohmann et al. (2013) describe the statistical analysis in more detail to evaluate linear temperature trends between 6 and 0 kyBP. Only those records are considered that have at least 10 incorporated values. The general temperature pattern is a warming in the tropics and cooling predominates in mid- and high-latitudes of the North Atlantic Ocean. In many regions, such as the North Atlantic Ocean, the Mediterranean Sea, the northern Indian Ocean, and the western North Pacific Ocean, there is a good agreement between the model and alkenone data with respect to the spatial pattern of the temperature trend (Fig. 8). The amplitudes of recorded and simulated temperature trends often differ, with proxies generally showing larger SST changes during the Holocene.

Fig. 8
figure 8

Holocene temperature trend (linear trend from 6 kyBP to pre-industrial) of the annual mean ECHO-G model output (see text for the details). The circles and squares localize the alkenone and Mg/Ca records, respectively. The colors that fill the circles/squares show the temperature trend they record. Units are ∘C per 6000 years

Full size image

4 Discussion

Our analysis indicates that on orbital time scale part of climate variability is related to the modulation of coherent patterns of climate variability. The effect of insolation forcing on atmospheric dynamics is isolated in the model set-up, and we have neglected variations in many subsystems of the Earth’s climate system relevant on these time scales, e.g., the ice sheets, carbon cycle, and vegetation. The results should be interpreted in a strict climatic sense only for interglacials. In the composite analyses (Fig. 3), we have taken the whole period from the model output in order to get a proper statistics and to better understand the associated mechanisms. The main finding is that the orbital forcing modulates the AO/NAO, NAM, and SAM modes affecting climate.

We focus on the boreal winter season since it strongly determines the long-term signal in SST. The persistence of SST anomalies is due to the fact that upper-ocean temperature anomalies created over a deep mixed layer in winter are preserved in the summer thermocline and reappear at the surface in the following winter. As discussed by (Joussaume and Braconnot 1997), a precise comparison of astronomical seasons under different orbital parameters would be to define the length of the months on the basis of their current angle on the Earth’s orbit. We find that the slightly different phasing of the winter season has no effect on our composite analyses.

The mechanism for the Northern Hemisphere dynamics is related to tropical convection. The convection anomaly and its associated diabatic heating directly drives the atmospheric circulation. The persistent tropical anomalies modulate the position of the storm track affecting the planetary-scale waves. Consistently, analysis of observational data indeed shows that the Pacific storm track is displaced poleward prior to the onset of the negative NAO phase (Franzke et al. 2004). In this context, it is worth mentioning that the ENSO-NAO relationship is non-stationary for the instrumental record (Rimbu et al. 2003b) and can have opposite trends for the last decades than for the entire observational period (Hoerling et al. 2001; Rimbu et al. 2003b). Hoerling et al. (2001) found a positive correlation of positive SST anomalies in the tropical Pacific (i.e., El Niño conditions) with the NAO when considering the last decades. However, negative SST anomalies in the tropical Pacific (i.e., La Niña conditions) are statistically significantly related to a SLP pattern in the North Atlantic similar to the positive phase of the NAO when the entire observational period is considered in the analysis (Pozo-Vázquez et al. 2001). Interestingly, the same AO/NAO-ENSO relationship holds for variations on orbital time scales.

It is very interesting to compare the ENSO modulation with other findings. Predominant La Niña conditions in the tropical Pacific during the early Holocene are suggested also by a simulation with an intermediate complexity ENSO model (Zebiak and Cane 1987) forced by variations in heating due to orbital variations in seasonal insolation. It shows weaker ENSO activity in the early-to-mid Holocene than in the late Holocene (Clement et al. 1999). Another model simulation shows that tropical Pacific temperature and circulation patterns during 6 kyBP are similar to those observed at the present-day La Niña period (Kitoh and Murakami 2002).

The astronomically driven AO/NAO in the model is in line with the analysis of the seasonal signal in an atmospheric circulation model (Hall et al. 2005). A sizable tropical SST response to the precessional component of Milankovitch variations is also seen in the three-dimensional global coupled ECBilt model (Tuenter et al. 2003), thus, indicating that the rectification effect seen by Clement et al. (1999) is not model dependent. It has been further argued that the precessional component of Milankovitch variations may be relevant also for deglaciation processes through low-to-high latitude teleconnections (Rodgers et al. 2003).

Interestingly, the SAM is in phase with the NAM. The SAM contributes a significant proportion of Southern Hemisphere mid-latitude circulation variability on many time scales (e.g., Hartmann and Lo 1998; Thompson and Wallace 2000). A negative SAM might have consequences for the large-scale ocean circulation and carbon cycle. A negative SAM is associated with an equatorward shift of the storm track and weaker westerlies above the Southern Ocean, less insulation of Antarctica and a warming of the Antarctica and the surrounding seas. These effects were documented in model experiments (Toggweiler and Samuels 1993; Sijp and England 2008; Wei et al. 2012). As a logical next step, future studies should include the interactions of atmospheric dynamics with other climate components such as the ocean carbon cycle, possibly amplifying the climate response as detected in our simulation.

Obliquity explains most of the variance in the annual insolation, and the effect is symmetric between the hemispheres but antisymmetric between the tropics and high latitudes. As a model for paleoclimate variations (Imbrie et al. 1993), quite often a seasonal index of insolation is taken (e.g., 65 ∘ N, 21st June), which provides already a strong linearity.

The precessional forcing (∼20 ky variations) largely dominates the seasonal effects, but is zero for local annual mean insolation forcing (Milankovitch 1941; Berger 1978; Berger et al. 1993). Short et al. (1991) used a linear 2D energy balance model to study the spatially resolved temperature response of the last 800 ky. Because of its linearity, the precessional response is absent in the annual mean response. A rectification can arise from many physical processes within the climate system, for example, a dependence of ice cover only on summer maximum insolation (e.g., Huybers and Wunsch 2003). A seasonal template model (Laepple and Lohmann 2009) provides a theoretical framework for the non-linearities indicating how the precessional frequencies are obtained. The precessional signal can be only due to non-linearities to the forcing (such as stronger response in summer than winter) or non-local responses (such as teleconnections as evaluated here). The strong boreal winter response to the tropical insolation forcing might be a key feature to multi-millennial climate variability and might even play a role for glacial-interglacial transitions. In our analysis of the glacial inception (120 to 115 kyBP), it seems that the boreal winter signal is at least in the same order of magnitude than the pure boreal summer signal. Both are acting in the same direction and provide synergistic memory effects: There is a delayed response of sea ice to summer insolation, there is a strong effect of winter SST on summer temperatures due to the ocean’s capacity of storing heat. Therefore, the atmospheric dynamics plays a key role on multi-millennial time scales. Future studies should include the interactions of atmospheric dynamics with other climate components such as land ice. As an example, the annual ablation and therefore the response of the glacial mass balance are sensitive to the integrated summer insolation (Huybers 2006; Huybers and Denton 2008).

For the mid-Holocene, high obliquity results in more high-latitude summer insolation at the expense of low-latitude summer insolation. The proxy-derived SST estimates indicate a similar pattern as the simulated annual mean SSTs supporting our finding of a multi-millennial AO/NAO and NAM variability. A necessary next step is also the comparison with Southern Hemisphere data. However, a systematic comparison with proxy data is beyond the scope of the paper. In several aspects, climate model and reconstructed temperature trends are to a large degree only qualitatively comparable, thus, providing a challenge for the interpretation of proxy data as well as the models sensitivity to orbital forcing (Lohmann et al. 2013). There is an additional difficulty when comparing the different frequencies in the model with data (e.g., Huybers and Wunsch 2003). The seasonality inherent in many climate proxies will produce precession-period variability in the records independent of any precession-period variability in the climate. In order to examine the precessional climate effects as evaluated here in the model, one shall have a deep insight into the recorder systems.

5 Conclusions

Applying an acceleration technique for Milankovitch forcing (Lorenz and Lohmann 2004), the modulation of large-scale climate modes on orbital time scales is addressed. Simulated Northern Hemisphere winter temperatures at mid-to-high latitudes are strongly linked to variations in the atmospheric dynamics on orbital time scales. We find that a the negative phase of AO/NAO (and negative phases of NAM as well as SAM) is associated with an increase in tropical convection. These atmospheric modes are modulated by tropical insolation during the boreal winter season associated to the Earth’s precession cycle. In contrast to boreal winter, the summer is largely affected by the local radiative forcing.

How can the variations of the climate modes on orbital timescales be understood? Given the relatively short time scale of the atmospheric variability, the atmospheric modes are modulated on time scales longer than that of the life cycles of the teleconnection patterns. The modulation is due to remotely forced stationary waves (Hoskins and Karoly 1981) which can be viewed as a boundary-value problem.

For the Eemian, the changed atmospheric circulation reflects a temperature pattern similar to the positive phase of the AO/NAO and mild winters in Europe (Zagwijn 1996; Aalbersberg and Litt 1996; Klotz et al. 1996). A large part of the strong winter cooling is due to the change in Northern Hemisphere circulation associated to a weakened Icelandic Low. It can be speculated that the glacial inception is due to the direct effect of insolation and atmospheric dynamics, i.e., summer cooling in combination with a negative phase of the AO/NAO during winter.

Coral data from the northern Red Sea suggest a modulation of AO/NAO for the Eemian and late Holocene climate (Felis et al. 2004). Our multi-proxy mapping effort reveals furthermore contrasting SST trends at one latitude. The spatial heterogeneity of the climate evolution at the end of the last interglacial and the Holocene trends supports a modulation of large-scale climate modes on orbital time scales.

Additionally to the dynamics caused by Earth’s orbital parameters, multi-centennial cooling events are reported for the central European continent based on an annually resolved, layer-counted record (Sirocko et al. 2005). More high-resolution proxy data are necessary in order to evaluate the spatial pattern related to such events. Future studies will also examine common and different mechanisms of climate variations during the Holocene and the last interglacial. Given the apparent involvement of the atmospheric dynamics on orbital time scales, mechanisms triggering the atmospheric modes should therefore not only provide interpretation of proxy data but also encompass and connect remote regions and identify their effects on long-term climate variability.