Antarctic Ice Core Harmonic Analysis

1. Introduction

The relationship between atmospheric CO₂ and temperature over glacial-interglacial cycles has been extensively studied using ice core records. A central question concerns causality: does CO₂ drive temperature through radiative forcing, or does temperature drive CO₂ through effects on ocean solubility, biological productivity, and terrestrial carbon storage?

Previous studies have established that CO₂ variations lag temperature changes during glacial terminations by several hundred to a few thousand years (Humlum et al., 2013; Grabyan, 2025). However, the precise phase relationship varies across different studies depending on methodology. Here we apply harmonic analysis with a novel period-selection algorithm to Antarctic ice core data spanning the last 350,000 years.

2. Data and Methods

2.1. Data Source

We analyze Antarctic ice core data converted to proxy temperature from the EPICA Dome C record (Jouzel et al., 2007), comprising 4,752 data points covering the interval from 350,245 years before present to 1911 CE. The corresponding CO₂ record (Bereiter et al., 2015) extends to 2001 CE. To avoid potential complications from apparent changes in interglacial maxima around 450,000 BCE (Higginbotham, 2025), the data interval used in the fit is from 350,000 BCE to present, or to 1850 CE for CO₂ to exclude anthropogenic influence. The Mid-Pleistocene Transition (MPT) associated with changes around 450,000 BCE is discussed below with respect to the temperature fit.

In the tables and figures that follow, negative dates will correspond to BCE while positive dates will correspond to CE.

Temperature: EPICA Dome C deuterium-derived temperature proxy (Jouzel et al., 2007), EDC3 chronology, 4,752 data points from 350,245 years BP to 1911 CE.

CO₂: Antarctic composite record (Bereiter et al., 2015), AICC2012 chronology. Data truncated at 1850 CE to exclude anthropogenic influence. The fit involved 1130 data points from 350,000 BCE to 1850 CE.

CH₄: EPICA Dome C methane record (Loulergue et al., 2008), EDC3 gas age chronology, matching the temperature chronology. The fit involved 1418 data points from 350,000 BCE to 1824 CE.

2.2. Harmonic Fitting Methodology

We fit each dataset to a harmonic model of the form:

X(t) = C + Σᵢ Aᵢ sin(2π(t − t₀)/Pᵢ + φᵢ)

where C is a constant offset, t₀ is the time origin (−175,000 years for all fits to allow direct comparison of phase), and each periodic component has amplitude Aᵢ, period Pᵢ, and phase φᵢ. For computational stability, the model is implemented using sine and cosine basis functions:

X(t) = C + Σᵢ [aᵢ sin(2π(t − t₀)/Pᵢ) + bᵢ cos(2π(t − t₀)/Pᵢ)]

where amplitude A = √(a² + b²) and phase φ = atan2(a, b). This parameterization avoids nonlinear optimization over phase.

2.3. Greedy Period Selection Algorithm

Rather than using predetermined periods or standard Fourier analysis, we employ a "greedy" algorithm: beginning with a constant term, iteratively add the period that most reduces residual variance, searching continuously over period space. The algorithm includes safeguards against ill-conditioning: minimum period separation, rejection of candidates causing excessive amplitude changes in existing periods, and amplitude ceilings.

2.4. Covariance Matrix Analysis and Period Refinement

To verify that selected periods form a well-conditioned basis, we computed the normalized covariance matrix for fitted parameters. Off-diagonal elements near ±1 indicate near-collinearity, which inflates parameter uncertainties.

Using temperature as an example

Initial analysis of the 59-period fit to temperature revealed problematic correlations:

• Periods 7,600 yr and 4,540 yr: r = −0.95
• Periods 8,900 yr and 9,000 yr: r = −0.70
• Periods 7,900 yr and 8,000 yr: r = −0.61
• Periods 23,000 yr and 22,150 yr: r = −0.51

We systematically removed the higher-index period from each problematic pair: 4,540, 9,000, 8,000, and 22,150 years. The final 55-period basis shows maximum off-diagonal correlations of |r| ≈ 0.15 for temperature and |r| ≈ 0.52 for CO₂—acceptable levels for stable parameter estimation.

Table 1. Fit quality before and after covariance refinement.

Configuration	R2 (Temperature)	Max |r|
Original (59 periods)	0.9539	0.95
Refined (55 periods)	0.9517	0.15

Table 1. Fit quality before and after covariance refinement.

Configuration	R2 (Temperature)	Max |r|
Original (59 periods)	0.9539	0.95
Refined (55 periods)	0.9517	0.15

The minimal R² change (0.0022) confirms these components were largely redundant.

Where CO₂ and CH₄ were independently fit, the same procedure was followed to obtain a good R² with acceptable |r|.

2.5. Uncertainty Estimation from Covariance Matrix

The unnormalized covariance matrix provides variances for each fitted parameter. Standard errors (SE) on sine and cosine coefficients propagate to amplitude and phase uncertainties via:

SE(A) = √[(a/A)² var(a) + (b/A)² var(b)]

SE(φ) = (180/π) × √[(b/A²)² var(a) + (a/A²)² var(b)]

Phase lag uncertainties combine temperature and CO₂ phase uncertainties in quadrature: SE(lag) = (P/360°) × √[SE(φ_T)² + SE(φ_CO₂)²].

2.6. CO₂ Fitting with Periods Found for Temperature

We fit CO₂ data to the same 55 periods with constant, allowing only amplitudes and phases to vary. Data after 1850 CE was excluded. The truncated fit achieves R² = 0.964, with maximum off-diagonal correlation |r| = 0.52. See Appendix A for further comments.The fit involved 1130 data points.

2.7. CO₂ Fitting with Periods Found Independently

The fit to CO₂ with periods found independently, with data after 1850 CE excluded, used a constant and 31 periods achieving an R² = 0.958, with maximum off-diagonal correlation |r| = 0.47. The fit involved 1130 data points.

2.8. CH₄ Fitting with Periods Found Independently

The fit to CH₄ with periods found independently used a constant and 55 periods achieving an R² = 0.873, with maximum off-diagonal correlation |r| = 0.40 .The fit involved 1418 data points.

3. Results

3.1. Temperature Fit Quality

The 55-period fit achieves R² = 0.952, explaining over 95% of variance in the 350,000-year record. The fit tracks glacial-interglacial cycles, the Eemian peak, the last deglaciation, the Younger Dryas, and millennial-scale Holocene variations including the Little Ice Age. See Appendix B for more details and fit amplitudes.

Figure 1. Fit to temperature proxy from date -350,000 to date 1850.

Figure 1. Fit to temperature proxy from date -350,000 to date 1850.

Figure 2. Fit to Temperature proxy from the prior interglacial warm period to date 1850.

Figure 2. Fit to Temperature proxy from the prior interglacial warm period to date 1850.

Figure 3. Fit to Temperature proxy from date -25,000 to date 1850.

Figure 3. Fit to Temperature proxy from date -25,000 to date 1850.

Figure 4. Fit to Temperature proxy from date -12,000 to date 1850.

Figure 4. Fit to Temperature proxy from date -12,000 to date 1850.

Figure 5. Fit of CO₂ (periods from temperature fit) in ppm from date -350,000 to 1850.

Figure 5. Fit of CO₂ (periods from temperature fit) in ppm from date -350,000 to 1850.

Figure 6. Fit of CO₂ (periods from temperature fit) in ppm from prior interglacial warm period to date 1850.

Figure 6. Fit of CO₂ (periods from temperature fit) in ppm from prior interglacial warm period to date 1850.

Figure 7. Fit of CO₂ (periods from temperature fit) in ppm from date -25,000 to date 1850.

Figure 7. Fit of CO₂ (periods from temperature fit) in ppm from date -25,000 to date 1850.

Figure 8. Fit of CO₂ (periods from temperature fit) in ppm from date -12,000 to date 1850. There is an obvious separation of the fit from the data for the last 10,000 years for this case that is not in the fit where periods are found independently. The data itself show larger excursions between dates 1500 and 1850 .

Figure 8. Fit of CO₂ (periods from temperature fit) in ppm from date -12,000 to date 1850. There is an obvious separation of the fit from the data for the last 10,000 years for this case that is not in the fit where periods are found independently. The data itself show larger excursions between dates 1500 and 1850 .

3.2. Recovery of Milankovitch Periods

The algorithm independently recovered periods matching canonical Milankovitch cycles (Milankovitch, 1941):

3.3. CO₂ Fit Using Periods from Temperature Fit and Phase Analysis

The CO₂ fit, using periods from the temperature fit, achieves R² = 0.964. For phase comparison, amplitudes were normalized to positive values (adding 180° to phase when amplitude was negative). Phase lag = (Δφ/360°) × P, with negative values indicating CO₂ lagging temperature. See Appendix C for fit amplitudes.

The fit of CO₂ data to the same periods was done: (1) to determine if the same periods would fit the CO₂ data, (2) to allow direct comparison of phase for each period. Additionally, it turned out that subsequently fitting the entire set of data from ~800,000 BCE to 1911 revealed the Mid-Pleistocene Transition which is discussed below.

A fit of CO₂ data using the greedy algorithm with pruning will be presented and discussed below.

3.4. Phase Lag Results with Uncertainties Where Both Fits Use the Same Periods

Table 2. Phase relationships for dominant periods (55-period fit).

Period (yr)	T Phase (°)	CO2 Phase (°)	Lag (yr)	SE(lag)	Cycle
113,600	233.6	221.5	−3,818	±10,400	Eccentricity
40,200	163.0	151.9	−1,239	±5,100	Obliquity
65,000	158.7	146.5	−2,203	±13,300	—
23,000	177.7	167.7	−639	±5,500	Precession
29,600	156.2	176.0	+1,628	±8,600	—
26,600	152.6	155.4	+207	±13,800	—

Table 2. Phase relationships for dominant periods (55-period fit).

Period (yr)	T Phase (°)	CO2 Phase (°)	Lag (yr)	SE(lag)	Cycle
113,600	233.6	221.5	−3,818	±10,400	Eccentricity
40,200	163.0	151.9	−1,239	±5,100	Obliquity
65,000	158.7	146.5	−2,203	±13,300	—
23,000	177.7	167.7	−639	±5,500	Precession
29,600	156.2	176.0	+1,628	±8,600	—
26,600	152.6	155.4	+207	±13,800	—

The large formal uncertainties reflect the challenge of precisely determining phases for individual sinusoidal components in a multi-period decomposition. However, the consistent pattern of negative lags (CO₂ lagging temperature) across the three major Milankovitch periods provides collective evidence stronger than any single period alone. The complete phase lag analysis for all 55 periods is provided in Appendix E.

The fit of CO₂ data using the greedy algorithm with pruning (31 periods), and a fit of the CH₄ data using the greedy algorithm with pruning are discussed below. The phase issue for these fits is addressed with cross-correlation.

3.5. The Eemian CO₂ Plateau

The Last Interglacial reveals profound asymmetry: from −130,000 to −117,000 years, temperature declined ~7°C while CO₂ remained at 275–280 ppm. This 13,000-year plateau indicates CO₂ absorption during cooling is much slower than release during warming.

3.6. Mid-Pleistocene Transition (MPT)

Consistent with the (MPT) (Clark et al., 2006; Berends et al., 2021), Figure 9 indicates that the periods that fit the data well from date -350,000 to 1911 are not appropriate for the temperature variations prior to date -400,000. The R squared value for this plot was 0.7.

The dominant ~100,000-year eccentricity cycle that characterizes the late Pleistocene is largely absent in the earlier record, which was dominated by the ~41,000-year obliquity cycle. The fit using periods determined on data from date -350,000 and later appears to be almost orthogonal to the data prior to date -400,000.

Figure 10 shows how a running computation of R squared produces a sharp drop for this fit at date -400,000 (400,000 BCE).

Figure 10. Here the average squared deviation and the variance are computed from date 1911 into the past as running quantities. This allows the computation of R squared for the fit from a given horizontal axis date forward in time to 1911.

Figure 10. Here the average squared deviation and the variance are computed from date 1911 into the past as running quantities. This allows the computation of R squared for the fit from a given horizontal axis date forward in time to 1911.

Figure 11. CO₂ fit with periods found using the greedy algorithm followed by pruning to minimumize off diagonal correlation.

Figure 11. CO₂ fit with periods found using the greedy algorithm followed by pruning to minimumize off diagonal correlation.

Figure 12. Zoom on previous plot of Fit to CO₂.

Figure 12. Zoom on previous plot of Fit to CO₂.

Figure 13. Zoom on previous plot of CO₂.

Figure 13. Zoom on previous plot of CO₂.

Figure 14. Zoom on previous plot. Compare to Figure 8.

Figure 14. Zoom on previous plot. Compare to Figure 8.

Figure 15. CH₄ fit with periods found by greedy algorithm followed by pruning to minimize off diagonal correlation.

Figure 15. CH₄ fit with periods found by greedy algorithm followed by pruning to minimize off diagonal correlation.

Figure 16. Zoom on previous Figure.

Figure 16. Zoom on previous Figure.

Figure 17. Zoom on previous Figure. Note the fit to the Younger Dryas notch at around 10,000 BCE.

Figure 17. Zoom on previous Figure. Note the fit to the Younger Dryas notch at around 10,000 BCE.

Figure 18. Zoom on previous Figure showing more clearly the fit to the Younger Dryas.

Figure 18. Zoom on previous Figure showing more clearly the fit to the Younger Dryas.

Figure 19. The times series for temperature and CO₂ for intervals, 350,000 BCE -1850 CE, 175,000 BCE – 1850 CE, 120,000 BCE – 1850 CE, and 60,000 BCE – 1850 CE, were cross correlated. The results are plotted here showing a lag in CO₂ that decreases as the input is restricted to more recent data until leading temperature for the case of 60,000 BCE – 1850 CE.

Figure 19. The times series for temperature and CO₂ for intervals, 350,000 BCE -1850 CE, 175,000 BCE – 1850 CE, 120,000 BCE – 1850 CE, and 60,000 BCE – 1850 CE, were cross correlated. The results are plotted here showing a lag in CO₂ that decreases as the input is restricted to more recent data until leading temperature for the case of 60,000 BCE – 1850 CE.

Figure 20. The times series for temperature and CH₄ for intervals, 350,000 BCE -1850 CE, 175,000 BCE – 1850 CE, 120,000 BCE – 1850 CE, and 60,000 BCE – 1850 CE, were cross correlated. All cases show a lag with respect to temperature.

Figure 20. The times series for temperature and CH₄ for intervals, 350,000 BCE -1850 CE, 175,000 BCE – 1850 CE, 120,000 BCE – 1850 CE, and 60,000 BCE – 1850 CE, were cross correlated. All cases show a lag with respect to temperature.

Figure 21. Here CO₂ is plotted against temperature. This plot shows curvature that tends to flatten as the temperature rises.

Figure 21. Here CO₂ is plotted against temperature. This plot shows curvature that tends to flatten as the temperature rises.

Figure 22. Here CH₄ is plotted against temperature showing a linear trend.

Figure 22. Here CH₄ is plotted against temperature showing a linear trend.

3.7. CO₂ Fit with Periods from Greedy Algorithm

The CO₂ fit found by the greedy algorithm followed by pruning to minimize off diagonal correlation, achieves R² = 0.958 using a constant and 31 periods. The maximum off diagonal correlation is 0.47.

3.8. CH₄ Fit with Periods from Greedy Algorithm

This independent recovery validates both the methodology and the orbital pacing hypothesis for glacial-interglacial cycles.

The greedy algorithm independently recovered orbital periods for all cases:

3.9. Cross Correlation

The fit to temperature was evaluated at 10 year intervals from 350,000 BCE to 1850 CE. Similarly the fit to CO₂ by the greedy algorithm with pruning (31 periods and a constant) was similarly evaluated. These time series were cross correlated and the results are shown in Figure 19 .

Cross-correlation was computed as a function of time lag τ, where positive τ represents CO₂/CH₄ shifted forward relative to temperature – a positive lag. A peak in correlation at negative lag indicates that the second variable (CO₂ or CH₄) follows temperature; a peak at positive lag would indicate the reverse. The correlation values shown are negative because the fits were evaluated relative to their respective means, and both gases increase with warming (positive correlation with temperature change from the mean).

Similarly the fit to CH₄ was evaluated and cross correlated with temperature. For all cases plotted, CH₄ shows a lag with respect to temperature.

The consistency of the lag for CH₄ indicates that the time series agree on registration of their respective quantity with time. This suggests that the changing lag associated with the CO₂ plot is either real or the time registration of the CO₂ quantity has an error.

3.10. Cross Plots

The time series used for the correlations of section 3.9 were used to generate cross plots. The cross plot for CO₂ indicates a curvature that tends toward flattening with higher temperature. This may have an influence on the asymmetry discussed in section 3.5. The CH₄ cross plot shows a linear trend. High temperatures are associated with the relatively narrow interglacial peaks. This limits the data available at high temperatures in these plots.

This saturation behavior of CO₂ is consistent with temperature-dependent ocean solubility, where the capacity for CO₂ release diminishes as the ocean warms.

4. Discussion

4.1. Temperature Drives CO₂ at Orbital Timescales

The consistent lag of CO₂ behind temperature at orbital periods, for fits using the same periods, is evidence that temperature drives CO₂, not the reverse. Both variables fit with R² > 0.95 using identical periods, with CO₂ systematically lagging. This supports: (1) orbital forcing modulates insolation; (2) temperature responds; (3) CO₂ responds to temperature through ocean degassing, circulation changes, and terrestrial carbon storage. These findings are consistent with previous phase analyses (Humlum et al., 2013; Grabyan, 2025).

The independent recovery of Milankovitch periods by the greedy algorithm, without prior specification, provides strong validation that these signals are genuine features of the data rather than artifacts of the fitting procedure. The consistency across temperature, CO₂, and CH₄ records further supports the orbital pacing hypothesis.

4.2. Mechanisms for Asymmetric Response

Atmospheric Transport Limitation: During warming, dissolved CO₂ is immediately available for release at the air-sea interface. During cooling, atmospheric CO₂ must find its way back to the ocean surface—a rate-limited process.

Plant Stress Feedback: As CO₂ declines, C3 vegetation becomes stressed near 200 ppm, reducing uptake and creating a floor below which CO₂ cannot easily fall.

Sea Ice Expansion: Growing sea ice caps the ocean surface, reducing gas exchange area precisely when cooling is most rapid.

4.3. The Modern CO₂ Anomaly

The fit extrapolated beyond the data shows gentle Holocene oscillation around 270–280 ppm. Actual CO₂ departing above 350 ppm after ~1850 CE lies far outside the orbital envelope. This analysis establishes the baseline from natural forcing; whatever drives modern CO₂, it is not the orbital mechanics governing the past 350,000 years.

5. Conclusions

1. Greedy harmonic analysis fits Antarctic ice core temperature, CO₂ and CH₄ over 350,000 years with R² > 0.95, independently recovering Milankovitch periods (Milankovitch, 1941). Covariance matrix analysis guided refinement of the temperature fit, in particular, from 59 to 55 periods, ensuring stable parameter estimates.The refinement from 59 to 55 periods eliminated all correlations exceeding |r| = 0.50, reducing the maximum from 0.95 to 0.15 for temperature.

2. CO₂ lags temperature by 600–4,000 years at orbital timescales, supporting temperature-driven carbon cycle dynamics. This conclusion is based on phase analysis shown in Table 2 section 3.4. However the formal standard error computation indicates large possible errors. This conclusion is based primarily on consistency for eccentricity and Obliquity. Cross- correlation of temperature and CO2 shows CO2 lagging temperature when applied to the full ~352,000 years of input data. However as data is clipped toward recent time the lag shifts to a lead for CO2. For CH4 the cross-correlation consistently indicates that CH4 lags temperature.

3. The Eemian shows CO₂ elevated for ~13,000 years while temperature fell ~7°C, indicating asymmetric absorption/release kinetics. This asymmetric behavior may be the reason direct comparison of phase is difficult.

4. The modern CO₂ spike is anomalous relative to the 350,000-year orbital pattern – obvious for all existing ice-core data.

5. The Younger Dryas temperature event is fit quite well as can be seen in Figure 4. The fit involves 111 fit parameters but this corresponds to roughly one parameter every 3000 years which seems quite reasonable for the observed accuracy of the fit, especially given the normalized covariance indication of little correlation between periods. The Younger Dryas event is even more prominent in the CH4 data and is harder to fit but it is reasonably well represented as shown in Figure 18.

6. The R squared analysis shown in Figure 10 appears to unambiguously find the Mid-Pleistocene Transition and clearly justifies limiting the fit to the interval chosen.