Millennial variations in atmospheric CO₂ during the early Holocene (11.7–7.4 ka)

We present a new high-resolution record of atmospheric CO₂ from the Siple Dome ice core, Antarctica, over the early Holocene (11.7–7.4 ka) that quantifies natural CO₂ variability on millennial timescales under interglacial climate conditions. Atmospheric CO₂ decreased by ∼10 ppm between 11.3 and 7.3 ka. The decrease was punctuated by local minima at 11.1, 10.1, 9.1, and 8.3 ka with an amplitude of 2–4 ppm. Although the explanations of carbon cycle mechanisms remain uncertain due to insufficient paleoclimate records and model simulations, these variations correlate with proxies for solar forcing and local climate in the southeast Atlantic polar front, eastern equatorial Pacific, and North Atlantic. Additional CO₂ measurements using better-quality ice cores and carbon cycle models are needed to confirm the observation.

1 Introduction

Future climate and ecosystem changes due to the continual increase in atmospheric carbon dioxide concentrations caused by human activities are inevitable (IPCC, 2013). Understanding the links between the carbon cycle and climate becomes important for accurate projection of future climate change. Atmospheric CO₂ is controlled by carbon exchange with ocean and land reservoirs, and increased CO₂ in the future and consequent changes in the earth system will in turn impact CO₂ levels via feedbacks (Friedlingstein et al., 2006). Due to the limited duration of direct measurements of atmospheric CO₂, which only started in 1957 (Keeling, 1960), our understanding of the carbon cycle dynamics is limited on longer timescales. Air bubbles occluded in Antarctic ice cores allow us to reconstruct ancient air and may help us better understand the mechanisms that control atmospheric CO₂ (Ahn and Brook, 2008, 2014; Bereiter et al., 2012; Higgins et al., 2015; Lüthi et al., 2008; Marcott et al., 2014; Nehrbass-Ahles et al., 2020; Petit et al., 1999).

Understanding the carbon cycle during interglacial periods is particularly useful because climate boundary conditions are similar to those of the near future. Previous work on late Holocene CO₂ records shows centennial CO₂ variability linked with climate, but the control mechanisms remain unclear, in part due to the potential mixture of natural and anthropogenic sources and sinks (Ahn et al., 2012; Bauska et al., 2015; Etheridge et al., 1996; Goosse, 2010; Indermühle et al., 1999; Rubino et al., 2013; Ruddiman, 2003, 2007). By contrast, CO₂ records for the early Holocene (11.7 to 7.3 ka) should reflect only natural CO₂ variability due to a smaller human population (Ruddiman, 2003).

The early Holocene (11.7 to 7.0 ka) is known as a relatively stable period in comparison with glacial periods. Several authors have linked centennial to millennial variability in the early Holocene to changes in solar forcing, including studies of the eastern equatorial Pacific (Marchitto et al., 2010), North Atlantic (Bond et al., 2001), and Southern Ocean (Nielsen et al., 2004), with responses in proxy records at ∼11.1, 10.1, 9.1, and 8.3 ka linked to solar variability (Bond et al., 2001; Marchitto et al., 2010). A weaker (stronger) solar activity has been linked with increased (decreased) ice-rafted debris in the North Atlantic (Bond cycle), dominant El Niño-like conditions (La Niña-like conditions) in the eastern equatorial Pacific, weaker (stronger) Asian monsoons, expansion (reduction) of sea ice in the Southern Ocean, and colder (warmer) sea surface temperature in the Southern Ocean (Bond et al., 2001; Marchitto et al., 2010; Nielsen et al., 2004; Reimer et al., 2004; Vonmoos et al., 2006). However, it is not clear what mechanisms are involved (Bond et al., 2001; Darby et al., 2012; Marchitto et al., 2010).

Atmospheric CO₂ on millennial timescales is mainly controlled by exchange with oceanic reservoirs and terrestrial carbon stocks. Existing atmospheric CO₂ records from EPICA Dome C (Dome C) show little variability of atmospheric CO₂ on millennial timescales from 10.9 to 7.3 ka (Monnin et al., 2001, 2004). However, high-frequency signals might be muted due to gas-trapping processes at this low-accumulation site (Spahni et al., 2003).

In this study, we measured 99 samples of atmospheric CO₂ with ages between 11.7 and 9.0 ka from the Siple Dome ice core. This new record complements the existing Siple Dome CO₂ record for 9.0–7.3 ka (Ahn et al., 2014). With this record, we investigate the relationship between atmospheric CO₂ and climate variations on centennial and millennial timescales. Siple Dome benefits from an accumulation rate 4.2 times higher than at Dome C and 1.8 times higher than at Taylor Dome (Table 1). A conservative estimate for the width of the gas age distribution in the Siple Dome record gives ∼42 years for the early Holocene (Ahn et al., 2014). Thus, the Siple Dome ice core allows high-temporal-resolution and higher-quality gas data with a more precise age scale and signals that are much less muted by the gas-trapping process. The temporal resolution on average during the early Holocene reaches ∼30 years as compared to ∼80 years in the Dome C record.

Table 1 Glaciological characteristics of Antarctic ice cores.

2 Methods

2.1 CO₂ measurements

A total of 247 individual ice samples from 99 depth intervals were measured by needle cracker dry extraction and gas chromatography methods at Seoul National University (SNU) (see Fig. S1 in the Supplement). We adopted the well-established measurement methods from Oregon State University (OSU) (Ahn et al., 2009) with minor modifications including sharpening of the tips of ice-crushing pins to increase the gas extraction efficiency and use of a newer model Agilent 7890 gas chromatograph (GC).

Briefly, ice samples were cut and trimmed carefully with a band saw in a −21 ^∘C walk-in freezer at SNU. All visible cracks were removed to eliminate potential CO₂ alteration by trapping modern air. An ice sample of ∼ 8–10 g was placed in a double-walled vacuum chamber maintained at about −35 ^∘C using cold ethanol circulation between the walls of the chamber while flowing ultra-pure N₂ gas (99.9999 %) into the chamber. The ice sample was crushed in the cooled chamber by 91 steel needles moving straight up to down using a linear motion (bellows) vacuum feedthrough. The liberated air from the ice was collected for 3 min in a sample tube in a cryogenic system maintained at 11 K. The CO₂ mixing ratio was determined by the Agilent 7890A GC equipped with a flame-ionization detector, using a Ni catalyst which converts CO₂ to CH₄ before measurement. Sample air was injected into a stainless steel sample loop, and the extracted air from each ice sample was analyzed twice. The GC system was calibrated daily with a standard air tank (293.25 ppm CO₂, WMOX2007 mole fraction scale, calibrated by US National Oceanic and Atmospheric Administration, Global Monitoring Division). To examine the linearity of the GC, ice samples from five different depth intervals (CO₂ concentrations of 239–251 ppm) were analyzed with two different air standards (188.9 and 293.3 ppm CO₂, respectively). The average difference in the results using the different standards was 0.4±0.9 ppm (1σ) (Table S1 in the Supplement).

2.2 Age scale of the Siple Dome ice core records

The Siple Dome samples are placed on the improved Siple Dome chronology developed by Yang et al. (2017), which is aligned with the Greenland Ice Core Chronology, 2005 (GICC05) using the synchronization of CH₄ and δ¹⁸O_atm time series. Abrupt CH₄ changes have been shown to be synchronous within about 50 years with abrupt climate changes in Greenland during the last glacial period (Baumgartner et al., 2014; Rosen et al., 2014). Using this principle, abrupt changes in the composite Siple Dome CH₄ data were aligned with abrupt changes in δ¹⁸O_ice from the NGRIP ice core (North Greenland Ice Core Project members, 2004; Rasmussen et al., 2006) at the 8.2 ka event and end of the Younger Dryas (Yang et al., 2017). For the time period of 11.64–8.10 ka, ages were updated from the original chronology of Severinghaus et al. (2009) by interpolating the age offsets at the tie points (Yang et al., 2017). For the time intervals outside of 11.64–8.10 ka, the age difference was set constant with the difference at the closest tie point. The modified gas ages are younger than the Severinghaus et al. (2009) ages by less than ∼110 years.

3 Results

3.1 The new high-resolution CO₂ record during the early Holocene

We obtained 99 data points that cover 622.14–539.06 m at SNU, corresponding to 11.7–9.0 ka (Fig. 1). To extend the record to 7.4 ka, we made a composite data set using a previous CO₂ record from the Siple Dome ice core covering 9.0–7.4 ka measured by the needle cracker system at OSU (Ahn et al., 2014) (Fig. 1). Between two and six replicates (2.6 and 2.4 on average for SNU and OSU data, respectively) from individual depth intervals were analyzed. The standard error of the mean of replicates from the same depth interval was 0.8 and 0.5 ppm on average for SNU and OSU data, respectively, ranging from 0.01 to 1.75 ppm. The sampling resolution is ∼30 years for 11.7–9.0 ka and ∼15 years for 9.0–7.3 ka.

Figure 1 High-resolution atmospheric CO₂ records obtained from Siple Dome ice core, Antarctica, during the early Holocene. (a) Pink and blue circles are Siple Dome ice core records obtained at Oregon State University (Ahn et al., 2014) and Seoul National University (this study), respectively. Lines represent 250-year running means, and dotted lines represent 2σ uncertainties calculated from Monte Carlo simulation. For the simulation, we produced 10 000 different sets of CO₂ concentrations which vary randomly with Gaussian propagation in their uncertainties. (b) The green line indicates 250-year running means of the original Siple Dome CO₂ data processed by high-pass filtering at $\mathrm{1}/\mathrm{1800}$ yr⁻¹. The blue line indicates 2σ uncertainties of the 250-year mean value and cannot be used to interpret variations on shorter timescales.

To make a composite record of atmospheric CO₂, we tested for bias between the two data sets. Siple Dome samples are from seven depth intervals between 538.55 and 490 m. A total of 16 samples were analyzed at both laboratories (Ahn et al., 2014). The SNU measurements were higher than the OSU measurements by 0.3±0.7 ppm (1σ) on average, indicating that the SNU and OSU results agree well (Table S2). The small offset of 0.3 ppm was added to OSU data before combining them with the SNU results.

3.2 Comparison with existing CO₂ records for the early Holocene

The new atmospheric CO₂ record from Siple Dome was compared to the existing CO₂ data from Dome C measured using the needle cracker at University of Bern (UB) (Monnin et al., 2001, 2004) and the existing CO₂ data from the WAIS Divide ice core measured by the needle cracker at OSU (Marcott et al., 2014) (Fig. 2a). On multi-millennial timescales, the baseline levels of the Siple Dome and WAIS Divide CO₂ records (Marcott et al., 2014) are higher than those from the Dome C (Flückiger et al., 2002; Monnin et al., 2004) record (Fig. 2a and c). The CO₂ offset between the Dome C and Siple Dome ice cores is 3–6 ppm (Fig. 2a and c).

Figure 2 (a) Atmospheric CO₂ records. Red dots: atmospheric CO₂ record from Dome C ice core. Red line: 250-year running means of atmospheric CO₂ record from Dome C ice core. Blue dots: atmospheric CO₂ record from Siple Dome ice core. Blue line: 250-year running means of atmospheric CO₂ record from Siple Dome ice core. Green dots: atmospheric CO₂ record from WAIS Divide ice core. Green line: 250-year running means of atmospheric CO₂ record from WAIS Divide ice core. (b) The blue line indicates 250-year running means of the original Siple Dome CO₂ data processed by high-pass filtering at $\mathrm{1}/\mathrm{1800}$ yr⁻¹. The green line indicates 250-year running means of the original WAIS Divide CO₂ data processed by high-pass filtering at $\mathrm{1}/\mathrm{1800}$ yr⁻¹. The red line indicates 250-year running means of the original WAIS Divide CO₂ data processed by high-pass filtering at $\mathrm{1}/\mathrm{1800}$ yr⁻¹. (c) CO₂ offset between Siple Dome CO₂ record and other published CO₂ records. Red line: CO₂ offset between Siple Dome CO₂ record and Dome C CO₂ record. Green line: CO₂ offset between Siple Dome CO₂ record and WAIS Divide CO₂ record.

The offset between Siple Dome CO₂ data in this study and other CO₂ data sets could be related to differences in the analytical methods used to make the measurements. To examine the inter-laboratory analytical offset, several Taylor Dome ice samples were analyzed at OSU (Ahn et al., 2014). The OSU results were higher than those at UB by 1.5 ppm on average. Taking the analytical offset between OSU and SNU of 0.3±0.7 ppm (1σ) into consideration, the 3–6 ppm CO₂ offset between the Siple Dome record (measured at OSU and SNU) and Dome C or Taylor Dome (measured at UB) cannot be entirely attributed to experimental offset.

To compare the new record to the existing records on millennial timescales, we calculate the Pearson correlation coefficient between Siple Dome CO₂ and existing CO₂ records. For this calculation, we use the Siple Dome and existing CO₂ record, which were smoothed and high-pass filtered at $\mathrm{1}/\mathrm{1800}$ yr⁻¹ (see Sect. 3.3 for detailed information). The offsets between existing CO₂ records and our data are also calculated (Fig. 2c). We use 250-year running means of CO₂ records for this calculation.

The correlation coefficient between Siple Dome CO₂ and WAIS Divide CO₂ during 11.45–9.02 ka is 0.02 (p=0.28) (Fig. 2b). The CO₂ offset between the WAIS Divide record and Siple Dome record is quite random (Fig. 2a and c) because of scattering in the WAIS Divide CO₂ record during the early Holocene period. The WAIS Divide CO₂ data during the early Holocene were reconstructed from the ice just below the bubble clathrate transition zone (BCTZ). Previous studies raised an issue about the possibility of high-frequency noise of the atmospheric CO₂ record in the ice just below the BCTZ (Lüthi et al., 2010; Shackleton et al., 2019). This phenomenon might be related to the gas fractionation effect because of clathrate layering during bubble–clathrate transformation. Gas content starts to be fractionated in the BCTZ because of the differential permeation of gas species when bubbles have transformed to clathrates. CO₂ concentration in the first layer of clathrates is more enriched with higher bubble-to-clathrate permeation rates. Below the BCTZ, gas content slowly homogenizes again through molecular diffusion (Bereiter et al., 2009), which can cause high-frequency noise to the ice below the BCTZ. Thus, the WAIS Divide CO₂ data are not sufficient to discuss millennial variabilities of the early Holocene.

The CO₂ record from Siple Dome is roughly correlated with the CO₂ record from Dome C during 11.45–7.45 ka (r=0.42, p<0.001). We observe the CO₂ offset of 3–8 ppm in the 250-year running means. The CO₂ offset between the Dome C record and Siple Dome record decreases continuously from 11.7 to 7 ka with small variations at around 9.3 and 8.3 ka (Fig. 2). The small variations in the Dome C CO₂ record (1.4 ppm, compared to 3.0 ppm for Siple Dome) can be explained by the lower sampling resolution (∼80 years for Dome C vs. ∼20 years for Siple Dome) and a stronger damping effect on CO₂ concentration change at Dome C due to the slower gas-trapping process at Dome C (Spahni et al., 2003).

The millennial CO₂ variations in the ice cores could be attributed to different degrees of in situ CO₂ production in ice. The in situ production of CO₂ is caused by carbonate–acid reactions (Anklin et al., 1997; Barnola et al., 1995; Delmas, 1993; Neftel et al., 1988; Smith et al., 1997a, b) and oxidation of organic acids (Tschumi and Stauffer, 2000). Although Antarctic ice cores have relatively low concentrations of carbonates and lower site temperatures compared to Greenlandic ice cores (Tschumi and Stauffer, 2000), it is estimated that the in situ production of CO₂ for Antarctic ice cores is smaller than 1.5 ppm (Bereiter et al., 2009). If the chemical alteration is the main cause of the millennial-scale CO₂ variations, we may expect to observe CO₂ age offsets among different cores because of dissimilar ice age–gas age differences. However, no available data set supports this possibility.

To further evaluate the in situ CO₂ production, we considered potential reactions. First, we compared the CO₂ with non-sea-salt Ca (nssCa) content in the ice to check the carbonate–acid reaction in the ice. The concentration of nssCa is mainly controlled by dust delivery, but it can also be produced partially by the carbonate–acid reaction in ice. Thus, we examined the concentration of nssCa ions in the Siple Dome and Dome C ice. The nssCa records do not correlate well with the filtered millennial CO₂ variations in both Siple Dome ($r=-\mathrm{0.33}$) and Dome C (r=0.15) records during the early Holocene (Figs. S2 and S3). In addition, the nssCa trends in Dome C and Siple Dome ice do not agree (Figs. S2 and S3), but millennial CO₂ variations do. Second, we checked the CO₂ production by oxidation of organic compounds (e.g., $\mathrm{2}{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}+\mathrm{HCHO}\to \mathrm{3}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}+{\mathrm{CO}}_{\mathrm{2}}$) in ice (Tschumi and Stauffer, 2000). The Dome C site is located further from the ocean than Siple Dome, and we therefore expect lower organic content in the Dome C ice. Concentrations of organic compounds at our sampling depths are not available. However, the concentration of oxidant H₂O₂ in the top 2.5–100 m in the Siple Dome core is below the detection limit of ∼0.02 µM (McConnell, 1997), although 0.02 µM H₂O₂ still has the potential to produce CO₂ and can increase the mixing ratio in bubbles by 5 ppm given sufficient supply of organic compounds (Ahn et al., 2004).

In summary, the existing Dome C CO₂ records covering the early Holocene share similar trends in the Siple Dome CO₂ record despite an offset in longer-term means of a few parts per million. We note that CO₂ offsets of several parts per million among different ice cores are common features in different time intervals such as the last millennium (Ahn et al., 2012; Monnin et al., 2004; Rubino et al., 2019; Siegenthaler et al., 2005) and Marine Isotope Stage 3 (Ahn and Brook, 2008; Bereiter et al., 2012), although they share the same trends of CO₂ change on multi-centennial to multi-millennial timescales. Thus, it is likely that the millennial CO₂ variations during the early Holocene in the Siple Dome and Dome C cores reflect atmospheric CO₂ changes.

3.3 Atmospheric CO₂ variations on the millennial timescale during the early Holocene

Figure 1 shows the CO₂ record from Siple Dome during the early Holocene. CO₂ increased by ∼8 ppm between 11.7 and 11.3 ka and then decreased by ∼10 ppm from 10.9 to 7.3 ka. The rapid CO₂ increase at 11.7–11.3 ka might be associated with abrupt warming in the North Atlantic and abrupt strengthening of Atlantic Meridional Overturning Circulation at the end of the last glacial termination (Marcott et al., 2014; Monnin et al., 2001). The long-term CO₂ trend is generally similar to that of the major water isotope (δD) variations in Antarctic ice cores reflecting Antarctic temperature variations (Fig. S4).

The Siple Dome CO₂ record shows millennial variability of ∼ 2–4 ppm with local minima at 11.1, 10.1, 9.1, and 8.3 ka (Fig. 1). These variations resemble variability in other paleoclimate records that has been linked to solar cycle variations on these timescales (Figs. 3 and S5).

Figure 3 Comparison of atmospheric CO2 with climatic proxy records over the early Holocene. The records were smoothed at ∼250 years and high-pass filtered at $\mathrm{1}/\mathrm{1800}$ yr⁻¹. (a) Atmospheric CO₂ record from Siple Dome (in this study). Dotted lines show 2σ uncertainties calculated from Monte Carlo simulation. (b) ¹⁴C production rate from IntCal04 Δ¹⁴C data (Marchitto et al., 2010; Reimer et al., 2004). (c) ¹⁰Be flux record from ice core on the GICC05 timescale (Finkel and Nishiizumi, 1997; Marchitto et al., 2010; Rasmussen et al., 2006; Vonmoos et al., 2006). (d) IRD stacked records from the North Atlantic regions on an untuned calibrated ¹⁴C age model (Bond et al., 2001; Marchitto et al., 2010). (e) North Greenland Ice Core Project (NGRIP) ice core isotope ratio on the GICC05 timescale (Rasmussen et al., 2006). (f) Sea surface temperature from the eastern equatorial Pacific indicating El Niño-like or La Niña-like conditions (Marchitto et al., 2010). The data were radiocarbon dated by accelerator mass spectrometry (AMS), which was recalibrated by the Marine09 calibration curve (Reimer et al., 2009). (g) Sea surface temperature from the Polar Front of the Southern Ocean on the chronology of Mortyn et al. (2003) (Nielsen et al., 2004). (h) Sea ice presence from the Polar Front of the Southern Ocean on the chronology of Mortyn et al. (2003) (Nielsen et al., 2004).

To examine the relationship between atmospheric CO₂ and the other paleoproxy data sets on millennial timescales, the Siple Dome CO₂ record was smoothed and high-pass filtered at $\mathrm{1}/\mathrm{1800}$ yr⁻¹ due to two necessities. First, it is likely that high-frequency variabilities of the atmospheric CO₂ record (decadal-scale variations and centennial-scale variations) are high-frequency noise of the atmospheric CO₂ record. Thus, we smoothed data sets to eliminate high-frequency variability. Before making a 250-year running mean, we made a 1-year interpolation, because sample spacing between data points covering the early Holocene is not constant. Second, to eliminate multi-millennial drift of the CO₂ record, the data were high-pass filtered at $\mathrm{1}/\mathrm{1800}$ yr⁻¹, following previous methods by Bond et al. (2001) and Marchitto et al. (2010). The proxy records were also processed in the same way as the CO₂ record to remove high-frequency variability and long-term draft.

We evaluated uncertainties of the smoothed and high-pass filtered CO₂ record using Monte Carlo simulation. Random sampling was made from a probability distribution for each measured value and its standard deviation. We repeated this series of simulations 10 000 times, which is shown as 2σ in Fig. 1 (see the Supplement for detailed information).

We calculated correlation coefficients between the filtered CO₂ and climate proxy series to understand their relationship with atmospheric CO₂ (Fig. 3; see the Supplement for methods). To calculate correlation coefficients between records, we selected data from 11.45 to 7.45 ka. Correlation coefficients, their significance, and maximum correlation lags are shown in Fig. 4 and Table 2. The CO₂ record from Siple Dome is anti-correlated with the stacked IRD record in the North Atlantic (Bond et al., 2001) ($r=-\mathrm{0.49}\pm \mathrm{0.1}$, CO₂ time lag of 120±155 years), sea surface temperature (SST) record in the eastern equatorial Pacific indicating El Niño-like or La Niña-like conditions ($r=-\mathrm{0.41}\pm \mathrm{0.13}$, CO₂ time lag of 50±219 years) (Marchitto et al., 2010), and sea ice in the Southern Ocean ($r=-\mathrm{0.35}\pm \mathrm{0.17}$, CO₂ time lag of 190±228 years) (Nielsen et al., 2004). On the other hand, the CO₂ record is positively correlated with summer SST (SSST) in the Southern Ocean ($r=\mathrm{0.35}\pm \mathrm{0.17}$, CO₂ time lag of 52±228 years) (Nielsen et al., 2004). The results may imply a tentative link between atmospheric CO₂ variations and climate change on millennial timescales. The time lags might be caused by age uncertainties of the proxy records and/or response time of atmospheric CO₂ to climate change (Bauska et al., 2015; Bereiter et al., 2012; Carvalhais et al., 2014).

Figure 4 Correlation coefficients between CO₂ and proxies with CO₂ time lag calculated from Monte Carlo simulation. Vertical lines in black indicate zero time lag. Vertical lines in blue indicate maximum correlation coefficients between CO₂ and proxies with CO₂ time lag. (a) ¹⁴C production rate and atmospheric CO₂. (b) ¹⁰Be flux and atmospheric CO₂. (c) SST in the eastern equatorial Pacific and atmospheric CO₂. (d) IRD from the North Atlantic and atmospheric CO₂. (e) SST in the eastern equatorial Pacific indicating El Niño-like or La Niña-like conditions and atmospheric CO₂. (f) SI in the eastern equatorial Pacific and atmospheric CO₂.

Table 2 Correlation between Siple Dome CO₂ record and climate proxy records. Column (a) shows correlation coefficients between CO₂ and proxies with CO₂ time lags. Column (b) shows correlation coefficients between CO₂ and proxies without CO₂ time lag. “With MC” shows mean values from the simulations taking age uncertainties into account. “Without MC” is the classic calculation of correlation, without taking age uncertainty into account. Significance of the lag correlations was assessed against 1000 repetitions of the lag correlation calculation using synthetic data stochastically generated to have the same red noise characteristics as the original series.

The anti-correlations we find are between the Siple Dome CO₂ record and the ¹⁴C production rate ($r=-\mathrm{0.49}\pm \mathrm{0.12}$, CO₂ time lag of $-\mathrm{20}\pm \mathrm{148}$ years) and ¹⁰Be flux ($r=-\mathrm{0.52}\pm \mathrm{0.08}$, CO₂ time lag of 110±63 years). This suggests that CO₂ and solar activity co-vary on millennial timescales (Fig. 4 and Table 2). These observations imply that atmospheric CO₂ variations might be influenced by climate change driven by solar activity on millennial timescales during the early Holocene (11.7–7.0 ka) (Fig. 4 and Table 2).

There are two outliers at ∼11.08 and 10.83 ka, which are far from the 250-year running mean (Fig. 1). Since the two outliers can enlarge the amplitude of actual CO₂ change, the data were processed except for the two values (Fig. S7). The Siple Dome CO₂ record, except for two data points at ∼11.08 and 10.83 ka, was smoothed and high-pass filtered at $\mathrm{1}/\mathrm{1800}$ yr⁻¹. With these processed data, we calculated correlation coefficients between the filtered CO₂ and climate proxy series again (Table S3). The relationship between CO₂ data except for two outliers at ∼11.08 and 10.83 ka and climate proxies is similar to the relationship between original CO₂ record and climate proxies, which shows that two outliers do not highly impact our interpretation.

4 Discussion

Possible carbon cycle control mechanisms in the early Holocene

Understanding a link between climate variations and solar activity on millennial timescales during the early Holocene is important to decipher carbon cycle mechanisms. However, the climate mechanisms have not yet been deciphered. A possible mechanism is that changes of solar activities may impact stratospheric ozone concentrations, which can change stratospheric and tropospheric circulation patterns (Meehl et al., 2009). Higher solar activity may enhance the precipitation in the Intertropical Convergence Zone (ITCZ) and South Pacific Convergence Zone (SPCZ) (Meehl et al., 2009; van Loon et al., 2007). Consequently, the intensified moisture at those areas would increase trade wind strength and upwelling in the eastern equatorial Pacific region. These conditions would lead to La Niña-like climate states on a millennial timescale (Marchitto et al., 2010). This change in the eastern equatorial Pacific might have affected the North Atlantic (Darby et al., 2012).

If the CO₂ variations we observe are affected by solar variabilities via climate, a number of mechanisms could be involved, including the terrestrial or marine carbon cycles, or both. We discuss three possibilities here. First, a close relationship between CO₂ and climate proxies in Antarctica (Jouzel et al., 2007) on multi-millennial timescales (Fig. S4) suggests that CO₂ variations on these timescales might be principally controlled by Southern Ocean processes. Atmospheric CO₂ can be controlled by temperature and salinity in the ocean (the solubility pump); solubility of CO₂ is greater in cooler and fresh surface waters (Broecker, 2002; Takahashi et al., 1993). The formation of deep water occurs in polar regions with high water density, where surface waters are cold; thus, the oceanic uptake of atmospheric CO₂ through this mechanism is stronger in polar regions (Sigman and Boyle, 2000). We observed a tentative link between atmospheric CO₂ and summer sea surface temperature (SSST) from the polar front region of the southeast Atlantic on millennial timescales (Nielsen et al., 2004), which implies that lower SSST in the Southern Ocean might have led to the reduction of atmospheric CO₂.

Increased sea ice extent might have blocked release of CO₂ from CO₂-rich deep water to the atmosphere, and therefore decreased atmospheric CO₂ concentration as previously suggested for glacial–interglacial CO₂ variations (Stephens and Keeling, 2000). Our Siple Dome CO₂ record is negatively correlated with the sea ice extent in the Southern Ocean, although the sea ice extent reconstruction shown in Fig. 3 represents only the east Atlantic region of the Southern Ocean.

Oceanic processes associated with El Niño-like and La Niña-like climate variation could also impact the carbon cycle. Marine sediment cores from the eastern equatorial Pacific show that solar activity proxies are well correlated with El Niño-like and La Niña-like climate variations in the eastern equatorial Pacific SST proxy record (Marchitto et al., 2010). The eastern equatorial Pacific is the region where CO₂-rich deep water upwells. Increased upwelling during La Niña-like conditions and resulting increased CO₂ outgassing have been suggested for the CO₂ increase during the last deglaciation (Kubota et al., 2014). Siple Dome CO₂ is anti-correlated with SST in the eastern equatorial Pacific on millennial timescales (Fig. 2), which may imply that La Niña-like climate can lead to higher CO₂ values.

Terrestrial carbon is involved with photosynthesis and respiration in plants, and with soil respiration (microbial and root respiration). Thus, terrestrial carbon is mostly controlled by temperature and precipitation (Davidson et al., 2000; Mielnick and Dugas, 2000). On multi-millennial timescales, when temperature in Greenland increases from 10.9 to 7.4 ka, atmospheric CO₂ decreases. Expansion of vegetation in the Northern Hemisphere may partially contribute to the decrease in atmospheric CO₂ (Indermühle et al., 1999).

A recent high-resolution study for the last 1200 years shows that centennial CO₂ variability was mainly controlled by terrestrial carbon, most likely in the high latitudes of the Northern Hemisphere (Bauska et al., 2015). The stacked IRD from the North Atlantic may be used for an indicator of cool conditions in the North Atlantic (Bond et al., 1992, 2001). The strong relationship between IRD and atmospheric CO₂ indicates that colder climate in the North Atlantic may lower atmospheric CO₂ by impacting terrestrial carbon stocks during the early Holocene.

δ¹⁸O_ice from the North Greenland Ice Core Project (NGRIP) ice core (Rasmussen et al., 2006) indicating temperature in Greenland also reveals millennial local minima at similar time intervals as those of CO₂ (∼11.4, 10.9, 10.2, 9.3, and 8.2 ka); however, atmospheric CO₂ and temperature in Greenland are mismatched at the earliest early Holocene and ∼8.2 ka. Thus, there is no significant linear relationship between CO₂ and temperature in Greenland on millennial timescales, and our calculation indicates that CO₂ leads temperature in Greenland on millennial timescales, though the correlation is still too small to assume any relationship ($r=\mathrm{0.21}\pm \mathrm{0.07}$, CO₂ time lag of $-\mathrm{130}\pm \mathrm{63}$ years).

Temperature in Greenland during the early Holocene might be partially influenced by the internal climate system and/or by low-latitude solar forcing indirectly. Two main cooling events in Greenland are recorded at ∼11.4 and ∼8.2 ka (Rasmussen et al., 2007). The well-known 8.2 ka cooling event is mainly influenced by the collapse of the Laurentide ice sheet (Merz et al., 2015) rather than by solar forcing; when temperature was colder in Greenland at ∼11.4 ka, solar forcing was higher, not reaching a minimum until ∼11.2 ka. It is also elusive whether solar forcing has an influence on climate in Greenland at ∼11.4 ka (Mekhaldi et al., 2020). In short, a linkage between atmospheric CO₂ and climate change during the early Holocene remains uncertain due to insufficient paleoclimate records and model simulations.

In this study, we observed that atmospheric CO₂ is highly anti-correlated with the ¹⁴C production rate and ¹⁰Be flux on millennial timescales with CO₂ time lag during the early Holocene (Fig. 3). The local minima of atmospheric CO₂ highly match with the local maxima of the ¹⁴C production rate and ¹⁰Be flux (minima in solar activity) at ∼11.1, 10.1, and 8.3 ka. The phenomena might be related to large variations in solar activity. However, the relationship between solar forcing and atmospheric CO₂ is different at ∼9.1 ka. The ¹⁴C production rate and ¹⁰Be flux are positively correlated with CO₂ at ∼9.1 ka on sub-millennial timescales, indicating that atmospheric CO₂ was at a local minimum at ∼9.1 ka when solar forcing was relatively high.

We also check the correlation of CO₂ with solar activity during the last 2000 years on centennial timescales (Fig. S8). A positive correlation between solar forcing and atmospheric CO₂ is observed during the Little Ice Age (LIA). There are two periods in which sunspots were exceedingly rare. During the Maunder sunspot minimum (1647–1715 CE), total solar irradiance (TSI) was reduced by 0.85±0.16 W m⁻². Atmospheric CO₂ records from Antarctic ice cores commonly show a decrease trend during this period (Ahn et al., 2012; Monnin et al., 2004; Siegenthaler et al., 2005; Rubino et al., 2019). During the Spörer Minimum (1450–1550 CE), TSI record during this period also shows a decrease trend. However, atmospheric CO₂ decrease is not significant in Law Dome and EPICA Dronning Maud Land (EDML) records (Monnin et al., 2004; Siegenthaler et al., 2005; Rubino et al., 2019), while the WAIS Divide ice record shows a decrease during this period (Ahn et al., 2012) (Fig. S8). However, atmospheric CO₂ decreases drastically at ∼1600 CE when TSI shows a local maximum, which is similar to the relationship between solar forcing and atmospheric CO₂ at ∼9.1 ka. To conclude, it is vague how solar forcing is related with atmospheric CO₂ variations on millennial timescales.

Comparing the early and late Holocene requires attention due to different boundary conditions during these two periods and anthropogenic CO₂ during the late Holocene (e.g., Ruddiman, 2003, 2007). Variations in solar forcing are large on a centennial timescale during the early Holocene. Thus, the solar output effect might be enhanced since the climate system has not responded linearly (Mohtadi et al., 2016). However, due to a decrease in summer insolation and the small variation in solar forcing during 7–1 ka (Berger, 1978), solar forcing might play a less important role during the late Holocene. Further studies are needed to understand the relationship between atmospheric CO₂ and solar forcing on shorter timescales during the early Holocene with more proxy records and numerical models.

5 Conclusion

In this study, we present a 30-year-resolution CO₂ record during the early Holocene. Our data show that millennial atmospheric CO₂ variability of 2–4 ppm correlates with several climate proxies such as IRD in the North Atlantic, sea ice extent in the Southern Ocean, and El Niño-like conditions in the eastern equatorial Pacific, all of which appear to coincidently occur with solar activity minima (Bond et al., 2001; Marchitto et al., 2010; Nielsen et al., 2004; Reimer et al., 2004; Vonmoos et al., 2006). The relationships with the proxies are consistent with changes in several different mechanisms that could impact atmospheric CO₂ on millennial timescales including changing CO₂ outgassing from the Southern Ocean and the eastern equatorial Pacific and changing terrestrial carbon storage in the Northern Hemisphere. Our new observations may improve our understanding of the relationship between interglacial climate and carbon cycles on millennial timescales in the absence of anthropogenic CO₂ perturbations. Further study should focus on clearly deciphering the millennial CO₂ control mechanisms with improved paleo proxy records and carbon cycle models.