On a New Theory of Climate Interference for MISs and Glacial Terminations from Antarctica Ice-core Records: Interference Model

1. Introduction

The Mid-Late Pleistocene glacial-interglacial cycles are particularly relevant for understanding the dynamics of the past climate system and assessing the Holocene interglacial period and current climate change. However, the fundamental characteristics of climatic variability remain unclear [1,2,3,4,5,6]. The state of current knowledge is well summarised in the IPCC [7] report cited in Ellis and Palmer [4], which affirms that “Orbital forcing is considered the pacemaker of transitions between glacials and interglacials (high confidence), although there is still no consensus on exactly how the different physical processes influenced by insolation changes interact to influence ice sheet volume.”

During the Pliocene-Pleistocene, Northern Hemisphere (NH) glaciation gradually developed in four subtrends [8,9] that progressively lowered the global average temperature and ice-volume growth [10,11,12,13,14,15,16,17,18], which was followed by a relative stasis of ice growth and temperature cooling in the last 600 kyr; this last period is termed the Mid-Brunhes Oscillation (MBO) [8,19]. The trend of climate cooling was triggered by non-orbital long-term changes in atmospheric greenhouse gases (GHG) composition [8,20,21,22,23]. By superimposing orbital forcings on this long-term trend of climate deterioration that caused relevant mean climate state changes [8,9], spectacular glacial-interglacial cycles in the main frequency bands of precession (~ 21-kyr), obliquity (~ 41 kyr), and short eccentricity (~ 100 kyr) were determined, the ‘orbital pacemaker’ of climate [24,25,26]. The Mid-Pleistocene Transition (MPT) from ~ 1.2 to 0.7 Ma marked a fundamental change in the climate system state from 41-kyr to ~100-kyr cycles and the beginning of the classic Mid-Late Pleistocene ‘Ice Ages’ [8,9,10,11,27,28]. In this paleoclimatological context, the Earth’s climate in the last 800 kyr (i.e., the post-MPT period including the MBO recovery) is crucial for understanding the wide alternation between cold glacial periods and warm interglacial periods characterised by ~ 100-kyr dominant climate variability (51.6% variance) [19].

The terms ‘glacial’ and ‘interglacial’ underline the changes in ice-volume and sea-level that are associated with the appearance and disappearance of large NH continental ice-sheets [3,5]. Contributions have been made in the field of glacial/interglacial characterisation using a modelling approach [29,30,31,32,33,34]. This top-down approach includes the mathematical representation of atmosphere/ocean/sea ice, land surface and vegetation, ice sheets, carbon cycle, and astronomical forcings to simulate climate signals. Other studies investigated the climate variability using different proxies by reviewing the common features and differences between them to assess the pattern of glacial/interglacial strength, termination duration and timing, insolation characteristics, and intra-MIS variability [2,3,5,6]. Understanding the causes of these oscillations, which are also related to the big unresolved question of the MPT origin, requires not only the replication of a ‘bulk’ glacial cycle and an understanding of how slightly different external forcing can lead to a wide range of responses [3], but also the replication of the internal structure of glacial cycles [9]. The understanding of the climate system is not satisfactory from a dynamic perspective of how the climate system evolves from glacial to interglacial states and vice versa [5]. The long ice-core records obtained by the European Project for Ice Coring in Antarctica have reinforced the view that each glacial/interglacial cycle has individual patterns in amplitude, length and profile, which remain unclear [1,2,3,5]. There is no simple astronomical factor that controls the amplitude patterns among interglacial periods, which seem to arise, at least in part, from the patterns of carbon dioxide (CO₂) [5]. The latter interpretation was supported and extended to glacial stages by Viaggi [19] on both CO₂ and CH₄ patterns. However, an observational structural approach to achieve internal organisation of glacial and interglacial stages and substages is lacking. The saw-tooth shape of the Mid-Late Pleistocene glacial cycles (slow cooling and rapid warming) reflects changes in the internal dynamics of the climate system, as asymmetry was not found in any orbital forcing [11,35,36]. These internal dynamics need to be clarified. Moreover, the true nature of glacial ‘termination’ (TER) is still an open question. Glacial TER is strictly defined as the midpoint of the rapid and wide transition between Mid-Late Pleistocene glacial and interglacial climate conditions, as recorded in marine stable oxygen isotopes [37], and marks the interglacial onset [5]. In contrast to TER, glacial inceptions mark the transition from a relevant interglacial to a glacial state by decreasing temperature, ice sheet build-up, and falling sea level. Understanding their origin is even more difficult.

Viaggi [8,9,19,38] emphasised both the need to study palaeoclimate proxies by variance-driven decomposition spectral methods, and to study these forcing-related components as ‘responses’ (primary forcing perturbation signal + cospectral feedback responses) of the climate system overcoming the nominal-centric paradigm. Indeed, climate responses resulted from an exponential feedback-driven (albedo mechanisms, carbon cycle, etc.) energy-transfer [16,39,40,41,42,43,44,45,46,47,48,49,50,51,52] from an initial, low-energy, and orbitally paced insolation change [24,25,26]. This new viewpoint involves an interesting perspective change in climate studies from being nominal-centric to response-centric (forcing-related), thus providing a conceptual approach showing how interfering signals may produce patterns/features that resemble what we observe with new interpretative insights into the climate system.

The present study is the first to examine the internal structures for all glacial and interglacial MISs and TERs over the last 800 kyr, and tries to answer some of the open questions concerning the Mid-Late Pleistocene climate variability based on this new observational perspective. Owing to the full-resolution variance-driven decomposition of EPICA records (δD, CO₂ and CH₄) [5,53,54,55] obtained using singular spectrum analysis (SSA) [19], it has become possible to provide a novel quantitative structural interpretation of glacial-interglacial cycles and glacial terminations. This bottom-up approach uses the response components of EPICA records (mid-term, orbitals, and suborbitals) [9,19] to reconstruct the envelope of the thermal response through a physical interference model, with the goal of contributing to the understanding of the intensity, amplitude, asymmetry, and internal variability of 73 MISs and seven TERs by means of a new theory of climate interference. This research is structured into two companion studies: the present study focuses on the interference model and its second part will elucidate the interference patterns (that is, signal configuration at a given time) (in progress).

2. Materials and Methods

Benthic δ¹⁸O measures changes in global ice-volume and deep water temperature, which are controlled by high-latitude surface temperature [11,56]. In the present work, the benthic LR04 δ¹⁸O record of Lisiecki and Raymo [57] filtered by SSA at 1 kyr [8] was used and resampled at 0.5 kyr via cubic spline interpolation during the last 800 kyr for uniformity with the EPICA stack record [19]. The LR04 δ¹⁸O benthic stack was obtained from 57 oceanic series globally distributed over a wide latitudinal range (Atlantic, Pacific, and Indian Oceans) [57]. LR04 is a global average signal from which regional variability was smoothed. The LR04 stack was orbitally tuned to an ice model driven by insolation that occurs on 21 June at 65° N in La93(1,1) orbital solutions, with a corrected sedimentation rate [57]. The insolation phase and amplitude errors between the La93 and La2004 [35] orbital solutions are very low and can be considered negligible for the last 800 kyr [8]. The EPICA stack is the mean signal among the highly covariant standardised δD, CO₂ and CH₄ SSA-filtered records [19] recalibrated by PIWG [5] based on the AICC2012 age model [58] with an average uncertainty of 1700 ± 995 yr BP [59]. This is because CO₂ and CH₄ records shown a general delayed pattern (cross-spectral analysis) with respect to the δD signal up to the half-precession component (namely total 91.2 % variance), suggesting a climate system affected by feedback processes that are forced by initial temperature change on the Earth’s surface [9,19]. The EPICA smoothed stack (stack_s) is the Savitzky-Golay filtered signal by 4-order least squares polynomial fitting across a moving window of 15. This compensatesfor the higher resolution of the original EPICA records. The EPICA stack_s better approximates the global benthic LR04 δ¹⁸O signal [9] (Figure 1), the latter being intrinsically related to the global/high-latitude atmospheric temperatures [11]. Despite potential age discrepancies, there is a striking similarity in the glacial-interglacial variability between the EPICA stack_s and LR04 δ¹⁸O records, even though amplitude offsets exist possibly resulting from the deep-temperature component of benthic δ¹⁸O [56,60,61] and/or regional differences. Indeed, the cross-correlation function (CCF) indicates a Pearson correlation at a 0-lag of 0.87 (Figure 1b), which was highly significant based on the correlation test (Figure 1c). This correlation was the highest among the single EPICA δD, CO₂ and CH₄ filtered records in the range of 0.80 to 0.83 [9], and suggests a relevant role of GHGs feedback-related signals in transforming the temperature proxy in the response stack. The CCF exhibited the highest coefficient (0.91) at ‒5 lag number, equivalent to a δ¹⁸O ‘bulk’ lags of ~ 2.5 kyr. The finding demonstrates that the pattern of global glacial–interglacial cycles also arises from patterns of Antarctic atmospheric CO₂ and CH₄. Comparisons of glacial-interglacial models suggest that Antarctic temperature changes represent global-scale temperature variations [62]. The modelling results of Yin and Berger [33] also indicate that the variations in the annual mean temperatures averaged over the Earth and over the southern high latitudes were controlled by GHGs, whereas insolation played a key role over the northern high latitudes. These differences are due to different geographical configurations, which lead to more climatic responses to seasonal insolation forcing over the Arctic, and more global and annual GHGs forcing over the Southern Ocean [33]. Hence, the EPICA stack can be considered a proxy of the global/NH atmospheric temperatures averaged by feedbacks-related CO₂ and CH₄ covariant changes [9]. The ten EPICA stack components obtained by SSA offer the most complete (up to the Nyquist’s limit) spectral decomposition available of this valuable record and were attributed to the putative forcings by means of spectral analysis [19], nominal cross-spectral analysis (orbital solutions, see Table 6 and Figure 24 in [9]) and other considerations for the lesser-known suborbital band [19], including cross-spectral analysis of the Hallstatt component with a total solar irradiance Holocene ¹⁰Be/¹⁴C proxy (see sect. 4. of Additional file 1 in [19]). Some millennial suborbital components have an unclear (9.7 kyr cycle stack) or unknown (3.7 kyr cycle stack) origin and require further investigation in these little-studied bands. However, they are structural components of the EPICA record that cannot be ignored even though they account for only 2.43 % in variance. Therefore, where possible, direct links between physical forcings and the structural components of EPICA were evaluated [9,19]. It follows that EPICA stack components (mid-term oscillation, obliquity modulation cycle, short eccentricity, obliquity, precession, half-precession, sun-9.7 kyr, sun-6.0 kyr, ‘unclear’-3.7 kyr and sun-2.5 kyr) [19] approximate the forcing-related thermal responses of the climate system. This offers the opportunity to study MISs structure on the EPICA stack record by graphic correlation with the LR04 δ¹⁸O curve (Figure 2).

To avoid confusion in the variety of MISs nomenclatures, we used the nomenclature of Railsback et al. [63], in addition to some informal modifications for the purposes of the present study. The peak age of 73 MISs for the last 800 kyr were collected by a peak detection algorithm from the SSA-filtered LR04 δ¹⁸O benthic stack [8,57], and reassigned to the EPICA stack by graphical (visual) correlation to overcome differences in age model and samples resolution between the original records (Figure 2). Graphic correlation (peak-to-peak within the MIS time-ranges of Railsback et al. [63]) requires judgment to assess which features correspond to one another and to distinguish noise from signal features. For these reasons, SSA filtering of both the LR04 [8] and EPICA stacked records [19] were used to identify δ¹⁸O MIS peaks and capture the age of the correlated EPICA peaks. The EPICA stack_s is used as graphic ‘reference’ for peaks identification to compensate for resolution differences between the original records. The red-labelled MIS in Figure 2 are substages recognisable in the EPICA stack but barely discernible on LR04 (MIS-1b/1c/11d/11e/14d).

Figure 2. Peak-to-peak correlation panel between (a) LR04 δ¹⁸O SSA-filtered record [8,57] and (b) EPICA stack/stack_s [19]. (c) Original MISs nomenclature from Railsback et al. [63]. Peak-to-peak correlation within the MIS time-ranges of Railsback et al. [63] (vertical black bars). MIS nomenclature slightly modified (colored labels) for the purposes of the present study. Red label: MIS recognisable in EPICA stack and barely visible on LR04 (1b, 1c, 11d, 11e, 14d). Purple label: new informal MIS recognisable both in EPICA stack and LR04 (13b, 13c). Blue label: renamed MIS (13d, 13e). TER: glacial termination. Vertical dashed lines: records mean. Figure 2c is from [63] with permission.

Figure 2. Peak-to-peak correlation panel between (a) LR04 δ¹⁸O SSA-filtered record [8,57] and (b) EPICA stack/stack_s [19]. (c) Original MISs nomenclature from Railsback et al. [63]. Peak-to-peak correlation within the MIS time-ranges of Railsback et al. [63] (vertical black bars). MIS nomenclature slightly modified (colored labels) for the purposes of the present study. Red label: MIS recognisable in EPICA stack and barely visible on LR04 (1b, 1c, 11d, 11e, 14d). Purple label: new informal MIS recognisable both in EPICA stack and LR04 (13b, 13c). Blue label: renamed MIS (13d, 13e). TER: glacial termination. Vertical dashed lines: records mean. Figure 2c is from [63] with permission.

Some of these are new informal substages (MIS-1b/1c). The long-lasting MIS-13a of Railsback et al. [63] was split into two substages (purple label) recognisable in both the EPICA stack and LR04 (MIS-13b/13c); consequently, MIS-13d and MIS-13e were renamed (blue label). Considering the not simple definition of some terms in Quaternary climatology (glacial/interglacial, stadial/interstadial), it is appropriate to specify the terminology used in the present work. According to the American Commission on Stratigraphic Nomenclature [64], glacials represent episodes during which extensive continental glaciers developed and interglacials are intervals during which the climate was incompatible with a wide extent of glaciers. Further qualitative low-order subdivisions (substages) include stadials and interstadials, reflecting secondary advances (colder substages) and recessions or standstills (warmer substages) of glaciers, respectively [5]. Thus, these terms are strictly related to the relative mean temperature, ice volume, and sea level as reflected in the historical concept of MIS [65]. However, classification of MISs in terms of glacial/interglacial stages or stadial/interstadial substages is sometimes difficult or lacking, and is rather subjective [5,63,66]. In an attempt to provide a more objective criterion, the first categorisation of MISs on a quantitative and structural basis is presented in this study, as it is fundamental to the purpose of the present research. MIS intensity (or strength) is defined as the position of its peak in the proxy scale, which is usually referred to as meaning relative ‘temperature/ice-volume’ (cold/glacial, cold-temperate, temperate, warm-temperate, warm/interglacial, etc.). The MIS amplitude is calculated on the EPICA stack as the mean depth of its peak within the lower and upper branches, and implies a hierarchical concept (stage/substage). The categorisation of MISs in terms of stage/substage was performed with respect to the average amplitude and intensity values of EPICA stack oscillations during the pre-/post-Mid-Brunhes Event (MBE, ~ 430 kyr) time intervals. When both the MIS amplitude and intensity were greater or less than the respective pre- or post-MBE averages, the oscillation was categorised at the corresponding stage rank (interglacial or glacial, respectively). In all other cases, MIS was considered a substage (stadial or interstadial). The position of the substage peak relative to the dominant wave of the EPICA short eccentricity component was used to define the glacial (negative pulse) or interglacial (positive wave) contexts of the substage. Being dominant (51.6% variance), this thermal condition determined the relevant climate background for defining the glacial/interglacial context of MIS substages. Thus, MISs are ultimately categorised into stages/substages according to their amplitude, intensity, and background context (glacial stage, glacial interstadial, glacial stadial, interglacial stage, interglacial interstadial, and interglacial stadial). The aforementioned structural criteria provide a practical categorisation of MISs for the purposes of the present study.

TERs are diachronous climatostratigraphic boundaries because of the potentially variable mixing time of the ocean, resulting in local temporal offsets of up to ~ 4 kyr between the MIS boundaries [5]. As transitions, the identification and dating of TERs may be more controversial than those of maxima or minima [63]. In the present study, the ages of TERs were obtained from [57], ratified by the International Commission on Stratigraphy [67], and recalibrated on the EPICA stack as the midpoint between the age peaks of neighbouring MISs (Table 1). For dating TER-1 and TER-5 in the EPICA stack, informal changes in the MISs nomenclature were not considered for consistency with the literature. The EPICA midpoint ages of the TERs were approximated according to data acquisition at 0.5-kyr resampled time-series [19].

The interference model was evaluated on ten EPICA stack components using interference indices (simple sum, INT-I, and variance-weighted version, INT-I_w) and multiple linear regression model to reconstruct the envelope of the EPICA thermal response. The best resulting interference signal (INT-I_w) is segmented into interference classes to identify the homogeneous region of interference intensity. A response index (RI) is calculated as the difference between the number of EPICA positive (stack > 0) and negative (stack ≤ 0) weighted components to further investigate the interference model. RI measures the polarity alignment according to the time of the ten thermal components, regardless of their intensity. For example, an RI of 4 indicates that seven thermal waves are positive (warming) and three are negative (cooling). Three different smoothed versions (RI_s) by a Savitzky-Golay filter based on least squares 4-order polynomial fitting across moving windows of 20, 100 and 300, respectively, were applied to examine the mid-term (20)/long-term (100 and 300) behaviour of the polarity among thermal responses. A variant of the RI calculated for the five suborbital components only (RI_sub) quantified the polarisation of these cycles. The concept of the time alignment of signals is sometimes expressed as a positive/negative response index (PRI/NRI), that is, the percentage of positive (warming)/negative (cooling) components out of the total. Finally, the 80 reference ages (73 MISs and seven TERs) calibrated on the EPICA stack (AICC2012 age model, [5,58]) were used to collect a selected dataset of the ten forcing-related EPICA stack components to investigate the interference model of MISs and TERs (abbreviated the MIS dataset) by the INT-I_w class.

3. Interference Model of EPICA Stack Components

3.1. Interference Indices

Interference is a linear physical phenomenon in which two or more wave signals of even different amplitudes, frequencies and polarities are superimposed without affecting each other to form a resultant complex wave of greater, lower or the same amplitude [69]. Interference is an energy redistribution process. Figure 3 depicts the temporal interaction among the EPICA variance-weighted (Figure 3a) and unweighted (Figure 3b) SSA components of different origin, frequencies and strength (% variance) compared with the LR04 δ¹⁸O and EPICA stack records. Positive values of the stack components indicate warming climate polarity and vice versa. The weighted interference index (INT-I_w) of EPICA (Figure 3a) reconstructs at a given time the envelope of the resulting thermal response as the vector sum of the ten variance-weighted zero-mean thermal stacked responses (R) as follows:

INT-I_w = (4.4 × R_MT) + (4 × R_OM) + (51.6 × R_SE) + (19 × R_OB) + (8.4 × R_PR) + (3.8 × R_HP) +
(1.8 × R_S9.7) + (2.6 × R_S6.0) + (0.63 × R_UC3.7) + (0.12 × R_S2.5)

where [19]:

R_MT: mid-term oscillation stack (4.4%);

R_OM: obliquity modulation cycle stack (4%);

R_SE: short eccentricity stack (51.6%);

R_OB: obliquity stack (19%);

R_PR: precession stack (8.4%);

R_HP: half-precession stack (3.8%);

R_S9.7: sun-9.7 kyr cycle stack (1.8%);

R_S6.0: sun-6.0 kyr cycle stack (2.6%);

R_UC3.7: ‘unclear’-3.7 kyr cycle stack (0.63%); and

R_S2.5: sun-2.5 kyr cycle stack (0.12%)

Figure 3. Correlation panel among the ten EPICA stack (a) variance-weighted and (b) unweighted components of different origin, frequencies and strength compared with EPICA stack/stack_s and INT-I records (left). Positive values of the stack components indicate warming climate polarity and vice versa. The EPICA INT-I_w (a) and INT-I (b) (blue lines) are the variance-weighted and unweighted vector sum, respectively, of the ten zero-mean stack components (light blue and black lines). The values of INT-I are multiplied by 10 to obtain a graphic scale comparable to that of INT-I_w. A simple (unweighted) interference model of the climate components fails to describe the envelope of the thermal response of EPICA stack.

Figure 3. Correlation panel among the ten EPICA stack (a) variance-weighted and (b) unweighted components of different origin, frequencies and strength compared with EPICA stack/stack_s and INT-I records (left). Positive values of the stack components indicate warming climate polarity and vice versa. The EPICA INT-I_w (a) and INT-I (b) (blue lines) are the variance-weighted and unweighted vector sum, respectively, of the ten zero-mean stack components (light blue and black lines). The values of INT-I are multiplied by 10 to obtain a graphic scale comparable to that of INT-I_w. A simple (unweighted) interference model of the climate components fails to describe the envelope of the thermal response of EPICA stack.

Instead, INT-I (Figure 3b) is an unweighted version of INT-I_w as follows:

INT-I = R_MT + R_OM + R_SE + R_OB + R_PR + R_HP + R_S9.7 + R_S6.0 + R_UC3.7 + R_S2.5

Figure 3 shows how a simple (unweighted) interference model of the climate components fails to describe the envelope of the EPICA stack. This is because the stacked components are calculated as weighted mean by variance of the standardised δD, CO₂ and CH₄ signals [19], and this neutralises the original variance content. Indeed, the INT-I record (Figure 3b) does not reproduce the overall shape of the EPICA stack curves, whereas the weighted version describes the latter more accurately (Figure 3a). This is also demonstrated by the cross plots in Figure 4, where the EPICA stack shows a robust linear correlation (R² = 0.87) with INT-I_w (Figure 4a) and a low correlation (R² = 0.61) with the unweighted version (Figure 4b). The standard error of the estimate is 0.34 for INT-I_w and almost double that of the unweighted version (0.60). In addition, the histograms in Figure 5 show that the weighted-interference model is a better approximation of the EPICA stack record. The distribution pattern of INT-I_w (Figure 5b) is similar to that of the EPICA stack (Fig 5a), whereas the unweighted version shows a different histogram (Figure 5c). Relevant features are the strong positive asymmetry (warming) and the ‘truncation’ of the negative tail (cooling) of both EPICA stack and INT-I_w, along with the multimodal distribution. This type of intensity asymmetry (smaller glacial intensity than interglacials) of the post-MPT cycles is different from the temporal asymmetry that is known as ‘saw-tooth’ shape and is clearly visible in the 0-mean records in Figure 2b. The intensity asymmetry also relates to the key question of why the pre-MBE interglacials appear to be cooler (i.e., less warm) than the post-MBE interglacials [33]. These asymmetries are discussed in two dedicated sections (3.4 and 3.5).

Since the INT-I_w is a smoothed version of the EPICA stack, a multiple linear regression exercise was conducted to obtain a more accurate reconstruction of the record and evaluate the role of the INT-I_w. As expected, the multiple linear regression model of the weighted components provides an excellent reconstruction of the EPICA stack record, as demonstrated by the robust goodness-of-fit statistics (R² = 0.97; Std. err. = 0.18; F = 4528; Sig. = 0.000) (Table 2). This is not surprising because the variables in the model are the structural components of the original variable. Here, the weighted interference vectors (variables, composed in the INT-I_w) are the kernel of the multiple regression equation, an outcome that is not entirely obvious. The low standard errors increase predominantly in the suborbital components in progressive order (i.e., increasing frequencies).

These analyses demonstrate how Antarctic stacked records during the last 800 kyr can be described by a variance-weighted interference model (INT-I_w) of independent, forcing-related thermal responses that interfere at a given time according to their relative differences in frequency, amplitude and polarity. INT-I_w essentially expresses the interference kernel of the Antarctic surface thermal forcing, and approximates the global/NH atmospheric temperatures. Interference is a linear physical phenomenon, which does not disagree with the nonlinear nature of climate records being nonlinear in the responses that interfere at a given time [8,9,19,38]. Climate signals paced by forcing are modified exponentially in their response components through the interplay of positive feedback mechanisms [16,40,42,44,46,48,49,50], which transfer much of the energy of the climate system [8,9,19,38]. The nonlinear phase locking, in which the external Milankovitch forcing affects the amplitude and period of climate oscillations to satisfy the nonlinear resonance condition [70], is likely the mechanism linking weak forcing perturbations with the induced feedback signals to produce exponential responses. Similarly, the parametric resonance model of Caccamo and Magazù [52] introduces a set of positive feedback loops in which the climate system response is non-linearly amplified with a net positive gain. This leads to a self-reinforcing feedback mechanism that amplifies the effects of the small disturbances connected with Milankovitch’s cycles including eccentricity [52,71]. Interference is a different phenomenon involving response components at a given time through a process of stochastic energy redistribution to produce the bulk climate wave.

Twelve INT-I_w clusters by equal-width segmentation (20 unit) have been named to identify homogeneous regions of interference intensity centred on a 0-mean class in the relative range ‒10 to 10 (Figure 6).

The terms ‘positive’ or ‘negative’ interference used in the segmentation class legend (Figure 6b) refer to the absolute intensity (from very low to very strong) of the resulting warming or cooling waves, respectively. The interference model allows recalibration of TER ages according to the structural notion of the ‘cancellation interference’ of thermal responses between glacial and interglacial peaks (INT-I_w near-0 value) (Table 3, Figure 7). Although the age difference between the EPICA criteria is small, for some TERs (TER-1/4/5/7) the midpoint interference intensity is significantly different from the neutrality (such as TER-2/3/6, which are coincident) and experiences an unexpected relevant positive bias, as can be seen from the values of the INT-I_w. TER-1 and TER-4 belong to the VLPI class, whereas TER-5 and TER-7 belong to the HPI and LPI classes, respectively (Figure 7). Therefore, to study the interference structure even of the recalibrated TERs, they have been added to the MIS dataset as additional cases with the ‘INT’ subscript (87 total cases).

3.2. Polarity Alignment

This section discusses the time alignment of the EPICA components using the response indices (Figure 8) and their frequency distributions (Figure 9). The RI shows a quasi-normal distribution, as evident by the symmetry of the histogram and linearity of the Q-Q normality plot. In particular, the RI centres on the null value, indicating the maximum occurrence probability of 50% warming and cooling alignments (RI = 0). The finding suggests that the short-term interference of the independent thermal components is a Gaussian stochastic process governed by chance. Instead, at the mid-term (RI_s 20)/long-term (RI_s 100 and 300) scales, the stochastic interactions increasingly deviate from normality and exhibit a multi-modal distribution that reveals patterns of consistency, suggesting a possible effect on the features of MISs (Figure 9). Long-term polarity indices are minimum (cooling alignment) at ~ 680–660 kyr (MISs-16) and maximum (warming alignment) at ~ 420–380 kyr (MIS-11c/d/e) (Figure 8). In contrast, MIS-5e shows a relatively positive RI_s of 100 but a negative RI_s of 300, indicating a long-term antipolar pattern. MISs-7 are associated with long-term strong negative polarity (RI_s 300), which is opposed to weaker maximum of the RI_s 100. Interestingly, MISs-1 are associated with long-term negative alignment. These examples suggest complex interaction patterns among the components of the climate system.

3.3. Constructive and Destructive Interference

The classification of interference into ‘constructive’ or ’destructive’ is based on the relative amplitudes of the weighted components compared to the INT-I_w when the resulting wave is larger or smaller than the larger (absolute) positive/negative component, respectively. This implies that constructive interference reinforces the wider positive (warming) or negative (cooling) thermal waves at a given time, whereas destructive interference weakens the resulting wave. Figure 10 depicts the interference ‘heat map’ of the MIS dataset (87 cases) sorted by both INT-I_w class and INT-I_w to highlight the interference model results. Figure 11 summarises the interference model. The EPICA stack vs. INT-I_w cross-plots with colour superposition of the INT-I_w class, symbol shape by MIS/TER (Figure 11a), and constructive/destructive interference (Figure 11b) confirm that INT-I_w is a good approximation of the EPICA stack record (R² = 0.81), as shown in Figure 4a for 1597 cases. At first glance, interglacial stages result from high to very strong (73%) and intermediate (27%) constructive interference of warming responses (from IPI to VSPI: MIS-1c, MIS-5e, MIS-7e, MIS-9e, MIS-11c/e, MIS-13a, MIS-15a/e, MIS-17c, MIS-19c) related to global δ¹⁸O depletion, and appear to be compactly clustered towards medium/very high values of positive interference (Figure 10). Instead, the glacial stages appear more dispersed along the spectrum of negative constructive interference, with the exception of MIS-20a, which falls unexpectedly into the VLPI class. In fact, 64% of glacial stages are associated with high/very high values of constructive interference of cooling responses (VHNI and HNI: MIS-2, MIS-4, MIS-6a, MIS-12a/c, MIS-14c, MIS-16a/c, MIS-18e) related to global δ¹⁸O enrichment. The other 36% are characterised by low negative constructive interference (LNI: MIS-8a, MIS-10a, MIS-14a; VLNI: MIS-18a) and the unexpectedly low positive destructive interference of MIS-20a (VLPI).

In fact, MIS-20a (VLPI) appears to be an anomaly of underestimated negative INT-I_w, owing to destructive interference with a large warming response of short eccentricity. This anomaly may only be apparent because the EPICA record does not reach MIS-20c which may have a more pronounced negative envelope of interference, similar to MIS-6c/8c/10c/14c with respect to the corresponding MIS-6a/8a/10a/14a (Figure 8). Indeed, this shape of increased interference (warming direction) in late glacial MISs occurs in the cases of MIS-6a/c, MIS-8a/c, MIS-10a/c, and MIS-14a/c, although it is not as pronounced, suggesting it is a structural feature of the thermal forcing related to TERs. This topic will be explored further in an upcoming study. In contrast, low-intensity MISs, including 90% of substages, fall in the low-interference regions (‒30, 30) where dominant destructive interference minimises the thermal response (HDI: MIS-5b, MIS-7d, MIS-9a, MIS-11b, MIS-14d, MIS-15b/d, MIS-17a/b, MIS-19b; VLNI: MIS-3c, MIS-9b, MIS-18a/b/c; VLPI: MIS-5a, MIS-11a, MIS-13d/e, MIS-15c, MIS-20a) (Figure 11b,c).

Concerning TERs, TER-2/2_INT, TER-3/3_INT and TER-6/6_INT share the same dating and fall consistently in the HDI region with destructive patterns and ‘cancellation’ interference, whereas TER-1/4/5/7 dated by midpoint appear away from the near-0 interference region and are characterised by an unexpected constructive interference and positive bias (VLPI: TER-1/4; LPI: TER-7; HPI: TER-5) (Figure 7 and Figure 10). This observation is likely an artefact caused by the inaccuracy of the midpoint dating criterion determined by the asymmetrical pattern of glacial cycles. In general, the absolute differences in dating terminations between the two EPICA criteria are not relevant as they are contained between 1.5–2 kyr, with the exception of TER-5 which goes up to 5.5 kyr. On a conceptual/structural basis the differences become important for TER-1/4/5/7. Thus, TERs are most appropriately dated and described by the structural criterion of null-interference and clearly fall in the HDI region with destructive patterns. While constructive and destructive interference showed a clear frequency distribution pattern within the interference classes (constructive dominant in high-intensity regions; destructive dominant in low-intensity regions) as logically to be expected (Figure 11c), the RI shows a more heterogeneous distribution, reflecting the Gaussian stochastic nature of the short-term interaction among climate components (Figure 11d). Indeed, contrary to what one would expect, in the VHNI region there are no MISs with minimum RI (up to –6 only, or NRI = 80%, MIS-2/8c), while those with RI = –4 (NRI = 70%, MIS-6c/12a/18e) are more frequent, with one case having an RI = 4 (NRI = 30% only, MIS-6b). The HDI class shows the widest distribution of RI values ranging from –6 to 6, whereas in the VHPI region and above, there are positive RI values from 6 (PRI = 80%, MIS-11e) to a maximum of 8 (PRI = 90%, MIS-5e/9e/11c/19c). Only one case in the SPI class had a RI of 0, equivalent to a PRI of 50% (MIS-11d). From Figure 11e-f, it is interesting to note that the millennial-scale suborbital components (half-precession, sun-related 9.7/6.0/2.5 kyr cycles, and unknown-3.7 kyr cycle) are significant in determining the rank of MISs. This is demonstrated by the RI_sub showing the greater positive polarisation (3, 5) of such cycles more frequently in interstadials than in interglacial stages. The polarisation is substantially opposite in stadials compared to glacial stages, where higher frequencies of negative suborbital alignments are observed (RI_sub = ‒3, ‒5).

In summary, the response-centric interpretation of the climate in terms of physical interference shows complex stochastic interactions of both intensity and polarity among ten strength-weighted thermal waves spanning the mid-term, orbital (including short eccentricity), and suborbital frequency bands, leading to constructive or destructive energy redistribution with both glacial and interglacial polarity at different ranks (stage, substage). This overcomes the fact that the features and timing of nominal forcings and climate responses may be significantly different because of nonlinear transformations by intermediate feedback mechanisms and age model discrepancies.

3.4. Intensity Asymmetry

Interestingly, the skewness does not characterise the cyclical stacked components of EPICA which are basically symmetrical as shown by histograms (Figure 12) and asymmetry statistics (Table 4). Here, D'Agostino's test evaluates the degree of deviation of the symmetry of sample distribution from that of the theoretical Gaussian distribution. The bootstrapping procedure was carried out for deriving robust estimates of standard errors to calculate D'Agostino's test. The test shows that thermal components have a not significant skewness, while INT-I_w and EPICA stack exhibit a pronounced positive skewness (Table 4, Figure 5). This suggests that the intensity asymmetry results from the interference process. In quantitative terms, the interference-induced asymmetry variation measured from the short eccentricity stack is 336 %, a remarkable result.

To further verify this result, an exercise was conducted to develop INT-I_w by progressively adding climatic components according to their relative strength in descending order and measuring the asymmetry statistics each time (Table 5). The progressive increase in both the positive skewness and maxima owing to the summation of new components until stabilisation after half-precession is a relevant feature. In particular, the suborbital components excluding half-precession do not modify skewness because of their low quantitative impact [19]. The tendency to reduce the minima indicates an overall damping effect of the negative interference, reflected in the truncation of the glacial tails of both the INT-I_w and EPICA stack histograms.

This analysis demonstrates the key role of interference in transforming the intrinsic features of the thermal responses (like symmetry) leading to the characteristic distribution pattern of the EPICA stack (positive skewness and negative tail ‘truncation’). Considering the main MISs, interference acts on the structural framework of the dominant short eccentricity response mainly in terms of bipolar intensity amplification and very high variability, with a bias in the warming direction (interglacial), especially after the MBE (Figure 13, Table 6 and Table 7). Indeed, by measuring on the main glacial/interglacial MISs the percentage effect of interference on the short eccentricity response, an amplification of 539 ± 400 % (mean and std. error) is obtained for interglacials against a value of 58 ± 28 % for glacials (Table 7). Really, of 31 major glacial/interglacial MISs, only three show a destructive interference pattern (9.7%) (MIS-7d/15b/20a). The pre-MBE interglacials exhibit a mean amplification of 117 ± 78 %, whereas the observed intensification for the post-MBE interglacials is 803 ± 646 %. The average strengthening of the pre-MBE glacials is 26 ± 34 % compared with 89 ± 43 % for the post-MBE glacials. The outstanding warming bias of the post-MBE interglacials (MIS-1a/1c/5e/7a/7c/7e/9e/11c) is due to the interference of warming responses of relevant mean intensity in both high-variance bands (obliquity and precession, with the exception of short eccentricity) and suborbital components. The average contribution of suborbital cycles to climate warming was more pronounced in the post-MBE interglacials (2.4) than in the pre-MBE interglacials (0.9), although weak overall. This led to the highest mean INT-I_w (84.5) with respect to the pre-MBE interglacials (78.2, MIS-13a/15a/15e/17c/19c). Surprisingly, the mean warming response of the short eccentricity was higher in the pre-MBE interglacials (51.9) than in the post-MBE interglacials (43.1). The amplification effect of the interference on major glacial MISs was more moderate than that on interglacial MISs. Post-MBE glacials (MIS-2/4/6a/6c/8a/8c/10a/10c) show a slightly colder average INT-I_w (‒68.2) than in the pre-MBE glacials (‒56.4) (MIS-12a/12c/14a/14c/16a/16c/18e/20a). Lower cooling of the pre-MBE glacial MISs occurs because all the high-variance bands (short eccentricity, obliquity, and precession) exhibit less pronounced average cooling than the post-MBE glacials. In particular, the average cooling of the short eccentricity in the pre-MBE glacial MISs is mitigated by the robust warming during MIS-20a. If we exclude MIS-20a from the average, the difference in glacial interference between pre- and post-MBE MISs is further reduced (INT-I_w = ‒66.8 vs. ‒68.2, respectively). Unlike interglacials, the cooling caused by suborbital cycles is slightly more pronounced during the pre-MBE time (‒0.8 vs. ‒0.6). Of special interest are the MIS-7d/15b stadials, whose antipolar patterns of destructive interference (short eccentricity warming vs. obliquity/precession cooling) dampen the envelope of climate responses into the HDI class. A destructive interference pattern similar to previous patterns but with more pronounced antipolar intensities is that of MIS-20a. Here, the cooling of both robust obliquity and very strong precession is dampened by a robust heat wave with short eccentricity.

Thus, the interference model explains the key questions of why the post-MBE interglacials appear warmer than the pre-MBE ones and why the post-MBE glacials appear slightly colder than the pre-MBE glacials. However, interference acts as a process of stochastic energy redistribution of thermal responses at a given time. This, in turn, is the effect of nonlinear energisation operated by little-known intermediate feedback mechanisms [4,8,9,16,19,38,24,40,42,44,46,48,49,50]. The climate system is dominated by positive feedbacks capable of making extremely strong bipolar amplification with responses of up to 5‒7 times the orbital forcings [8,9,38]. This may explain the energy excess ‘paradox’ of the astronomically paced signals compared to the small energy of the orbital forcing, especially with respect to the eccentricity bands [10,26,31,52,72,73]. The energy paradox of Milankovitch's theory between the power of high-latitude NH summer insolation spectra and the variance of climate proxies can be understood if we change the climate interpretation from the nominal-centric to a response-centric perspective, with the feedback mechanisms in between [9]. The latter mechanisms make the difference because precession, although dominant in the NH summer insolation spectra, is a short cycle in which feedbacks have less time to transfer additional energy. This outstanding period-dependent connection is derived from the Plio-Pleistocene exponential relationships between the variance of orbital LR04 δ¹⁸O components vs. the mean δ¹⁸O exponential trend by time segments binned at 532 kyr [8]. Indeed, the exponential coefficient e^kx (x = time) for the δ¹⁸O semi-precessional component is 1.05 and increases with the orbital period; it is 1.78 for the precessional component, 1.88 for the obliquity component and 2.93 for the short eccentricity component [8]. These data indicate an increase in δ¹⁸O exponential growth from semi-precession to ~ 100-kyr eccentricity responses towards the Mid-Late Pleistocene as function of the long-term mean climate state (boundary condition) likely affecting the amplification effect of covariant feedback mechanisms linked to the forcing period [8]. Obliquity, although weaker in NH summer insolation with respect to precession, has a longer cycle that favours a more consistent (nonlinear) feedbacks-mediated thermal response. Therefore, an ice sheet that manages to survive the ablation of precession and obliquity peaks could be sensitive to feedbacks mediated by short eccentricity during the Earth’s long-term ice state (obliquity damping hypothesis) [9].

3.5. Time Asymmetry

The saw-tooth shape of Mid–Late Pleistocene glacial cycles is a well-known but poorly understood feature. This section investigates whether the physical interference of thermal responses plays a role in influenced this type of asymmetry during the last 800 kyr. The analysis was conducted by examining the cooling and warming durations of eight short eccentricity cycles using the EPICA stack weighed component (Figure 13). Cycles were defined between the EPICA stack short eccentricity maxima (from 1^st to 8^th in order of increasing age) and the La2010 eccentricity nominal solution [36] by capturing the ages of the maxima and minima with a peak detection algorithm (Table 8).

The maximum of the 1^st cycle of the EPICA short eccentricity stack, not yet reached at the present time, was estimated by an autoregressive linear forecast algorithm (40-order singular value decomposition) at ‒4.5 kyr in the future time. These data were compared to those obtained from structurally equivalent cycles of INT-I_w. Duration statistics were calculated for eight cycles for each record. The first unexpected result of this analysis is that the nominal short eccentricity cycles exhibit a slight mean time asymmetry, which is not reported in the literature [11,35]. Indeed, the nominal short eccentricity cycles show that the average duration of the cooling hemicycle (decreasing eccentricity) is 50.9 ± 14 kyr, while the warming hemicycle (increasing eccentricity) has an average time of 45.1 ± 11.2 kyr, resulting a mean time asymmetry (prolonged cooling) of 5.8 kyr. The EPICA short eccentricity stack exhibits an average cooling duration of 54 ± 14 kyr and an average warming time of 45.4 ± 7.1 kyr, with a mean time asymmetry of 8.6 kyr. Thus, the asymmetry in time of the EPICA short eccentricity response was amplified with respect to the nominal forcing by the 48% prolonged cooling. Moreover, the effect of interference measured on the INT-I_w is to lengthen the average cooling times to 59.1 ± 19.5 kyr and to shorten the warming times to 38.0 ± 15.7 kyr, thus bringing the time difference to 21.1 kyr prolonged cooling. Notably, these data show an intrinsic temporal asymmetry in the nominal short eccentricity and its EPICA response, as well as a strong interference-induced time asymmetry that prolonged the cooling vs. short eccentricity response by 145%. The 48% increase in the duration of the cooling hemicycle of the thermal response paced by short eccentricity forcing can be attributed to feedback mechanisms induced by an already asymmetric nominal hemicycle. Indeed, a longer forcing cooling hemicycle will give rise to a more prolonged feedback mechanisms compared to a shorter warming hemicycle, accentuating the intrinsic temporal asymmetry of the resulting thermal response. Thus, the overall effect of the increased temporal asymmetry (from nominal eccentricity = 5.8 kyr to INT-I_w = 21.1 kyr, i.e. 264%) which resulted in the saw-tooth shape of glacial cycles, is due to the joint action of intermediate feedback mechanisms with a strong contribution by physical interference. Finally, the average period of the short eccentricity cycle among the nominal forcing (96.0 kyr), EPICA response (99.4 kyr), and interference envelope (97.1 kyr) did not change significantly through feedback processes or the interference mechanism.

4. Relationships Among Orbital Forcings, EPICA Responses and Interference

As previously stated, the climate system consists of thermal responses paced by forcings and nonlinear changes through the synergistic action of internal feedback mechanisms which are both climate state-related (initial conditions) and scale-dependent. In other words, the intensity of the thermal responses is not strictly related to the intensity of the insolation/forcing but rather to the time (period) of the orbital disturbance affecting the covariant feedbacks [8]. This was probably because the determination coefficient (R²) between the nominal forcings and the related EPICA thermal responses was low (Figure 14), whereas the spectral cross-coherency relationships were more stable (Table 9). In fact, the cross-coherency was very high (eccentricity/short eccentricity = 0.94; obliquity = 0.98) or fair (precession = 0.67), whereas R² was low (eccentricity = 0.16; short eccentricity = 0.46; obliquity = 0.50; precession = 0.41). Note the very high coherency of the eccentricity/short eccentricity and the strong difference in R² between the nominal eccentricity (0.16) and its 400-kyr filtered-out version (0.46).

Assuming a certain margin of error due to the different age models, the phase delay of the EPICA obliquity and precession stacks with respect to nominal forcings (4.2-kyr and ~4.0-kyr, respectively) is in agreement with the putative cause-effect relationship and supports the action of feedback mechanisms that are lagged by definition [8,9,19]. On the other hand, the EPICA short eccentricity stack shows an advance of 4.3‒4.5-kyr, which seems uncredible. Instead, the phase relationships of the short eccentricity within the EPICA records (i.e. with the same age model) are consistent with an inertial climate system affected by feedbacks mechanisms since CO₂ and CH₄ lags δD [19]. The EPICA ‘anomalous’ result with the eccentricity nominal solution could be an artefact owing to the averaging of the lead-lag phase variability over a small number of short eccentricity peaks (16) during the last 800 kyr, perhaps in relation to possible inaccuracies in the age model and error propagation on long periodicity. Despite this inconsistency in the phase shift of the eccentricity nominal solution, and considering the very high cross-coherency, these data suggest that nominal short eccentricity played a key role in pacing the climate response over the last 800 kyr [8,9,19,24,40,49,74,75]. This is also supported by Figure 14, which shows the distribution of MIS interference classes which appear widely dispersed, especially for obliquity and precession (Figure 14c,d), but with some regularities for the short eccentricity response because of its relevant quantitative impact (MISs of the high positive interference class clustered towards high short eccentricity and viveversa; Figure 14b). On the other hand, the weaker short eccentricity external forcing should not cause any problem if the ~100-kyr response is driven primarily by internal feedbacks which transfer most of the climate system energy [8,9,19,24,40,49,52,74,75]. The absence of a univocal relationship between the strength of the astronomical forcing and the intensity of interglacials [5] can be explained and generalised to all MISs by the interference model, which shows that nominal forcing of the same strength corresponds to MISs of even very different interference classes resulting from the stochastic redistribution of the climate system energy through the physical interference of distinct thermal waves (Figure 14). The EPICA interference model includes the most complete suite of climate forcings in play (mid-term, orbitals, and suborbitals, including sun-related cycles), their nonlinear transformations through feedback mechanisms (climate responses), and their envelope by physical interference over time. If we change the interpretative paradigm from a nominal-centric view to the related climate responses [9], then the interference model may explain many issues of the Mid-Late Pleistocene climate variability. Climate interference is a cross-cutting theory that represents a quantitative model linking the non-linear responses of climatic components and their stochastic envelopes at a given time.

5. Conclusions

The main conclusions regarding the interference model of the EPICA thermal components are summarised as follows.

1)

The Antarctic stacked records during the last 800 kyr, which approximate global/NH atmospheric temperatures, are described by a variance-weighted interference model of independent thermal waves of different origins (mid-term oscillation, orbitals, and suborbitals) that stochastically interfere at a given time according to their relative differences in frequency, amplitude, and polarity.

2)

Climate interference over a short-term timescale is a Gaussian stochastic process governed by chance. At the mid-/long-term scale, the stochastic interaction of the responses increasingly deviates from normality and exhibits a multi-modal distribution that reveals patterns of time alignments, suggesting a complex interaction among the components of the climate system.

3)

Interglacial stages resulted from high to very strong (73%) and intermediate (27%) constructive interference of warming responses, and glacial stages derived by high/very high (64%) and low (36%) constructive interference of cooling responses, which are related to global δ¹⁸O depletion or enrichment, respectively.

4)

MIS substages are the result of positive (interstadial) or negative (stadial) polarisation of millennial suborbital thermal waves.

5)

Low-intensity MISs, including 90% of substages, fall in low-interference regions where dominant destructive interference minimises the thermal envelope.

6)

The interference model allows the recalibration of TER ages according to the structural notion of ‘cancellation interference’ of the thermal envelope. For TER-1/4/5/7 dated by the midpoint criterion, unexpected constructive interference leads significantly away from the neutrality region with a relevant positive bias, which is an artefact due to the asymmetry of glacial cycles. Therefore, terminations are most appropriately dated and described by the structural criterion of null-interference which is associated with destructive patterns such as TER-2/3/6.

7)

The intensity asymmetry of the EPICA stack record results from the interference process which, in turn, amplifies the bipolar envelope of nonlinear, symmetric thermal responses that interfere at a given time with a strong bias in the warming direction (interglacial, 539 ± 400 %), especially after the MBE (803 ± 646 %). This explains the key question of why the post-MBE interglacials appear warmer than the pre-MBE interglacials. However, even the post-MBE glacials appear slightly colder (89 ± 43 %) than the pre-MBE glacials (26 ± 34 %), albeit to a lesser magnitude.

8)

The duration analysis of short eccentricity hemicycles exhibits an unexpected, intrinsic mean time asymmetry (prolonged cooling) in the nominal short eccentricity (5.8 kyr) and its EPICA response (8.6 kyr, 48% variation), as well as an interference-induced asymmetry (21.1 kyr) that amplifies the cooling duration of short eccentricity response by 145%. Thus, the overall effect of the increased temporal cooling resulting in the saw-tooth shape of glacial cycles (from 5.8 kyr to 21.1 kyr, i.e. 264%) is related to the joint action of phase-locked intermediate feedback mechanisms, and especially by the strong amplification of the physical interference process.

9)

Nominal forcing of the same strength corresponds to MISs (interglacials and glacials) of even very different interference classes resulted from the stochastic redistribution of the climate system energy through the physical interference of feedbacks-amplified thermal waves. This explains why there is no simple relationship between orbital forcing and MISs intensity.

10)

The nominal short eccentricity played a key role in pacing the feedback-amplified climate response and its central action in the interference process over the last 800 kyr.