Northern Hemisphere atmospheric history of carbon monoxide since preindustrial times reconstructed from multiple Greenland ice cores

Hemisphere

. Locations, site characteristics and other relevant information for ice cores featured in this study. improve signal resolution. This section reviews the operation and recent improvements of CFA-based CO measurements. More details can also be found in Supplementary Information (SI :::: Sect. : 1).

System operation 120
The ice cores listed in Table 1 were analysed using continuous ice core melter systems coupled with online gas measurements.
The general principles of this analytical approach are provided by (Stowasser et al., 2012). Briefly, ice core sticks are cut at a 34 mm x 34 mm cross section and processed on a melterhead located in a cold room. The melterhead is composed of inner and outer collection areas with the inner area dedicated to sample collection. To prevent contamination, an overflow from the inner to the outer melterhead areas of >10% is created by lowering the sample pumping speed. The water and gas bubble mixture 125 is continuously pumped via a debubbler into a temperature-controlled gas extraction unit maintained at 30°C. The gas/water volume ratio of the sample is about 10% before the debubbler, and 50% after the debubbler. A fully degassed fraction of the melted sample is thus also available at the debubbler for complementary chemical analyses in the liquid phase. The gas is extracted along the sample line after the debubbler by applying a pressure gradient across a gas-permeable membrane, either Transfer-Line (Idex) or Micromodule 0.5×1 (Membrana GmbH) degassers. Then, the gas is dried by a custom-made Nafion 130 (Perma Pure) dryer. Finally, CO (and/or CH 4 ) mixing ratios are continuously monitored along the gas sample flow by a laser spectrometer.
Our samples were all analysed at the Desert Research Institute ::::: (DRI), Reno (NV, USA) during two analytical campaigns (Tunu13, D4, NEEM-SC, and NGRIP in 2013;PLACE in 2015). The NEEM-2011-S1 core was analysed at DRI in 2011, as described by Faïn et al. (2014). The PLACE core was also analysed at IGE ( ::: the ::::::: Institute :: of :::::::::::: Environmental ::::::::::: Geosciences ::::: (IGE, analysed by the laser spectrometer. However, it is not completely identical as it includes more components such as an additional peristaltic pump or extra lines. This pathway for synthetic standard analysis will be referred to as "calibration loop" here, but has been designated as "full loop" in previous studies (e.g., Rhodes et al., 2013;Faïn et al., 2014). Note that this "calibration 155 loop" should not be confused with the "Internal Loop" (Rhodes et al., 2013) a system to isolate the gas extraction system from the remainder of the melting system, while keeping a continuous gas extraction (see Fig. 1 from Rhodes et al., 2013).

CO blank of CFA-based analyses
The CO procedural blank was evaluated for each analytical campaign. Two different evaluation approaches were applied, which are described in SI :::: Sect. ::: 1.5. A relatively high CO blank of 35±7 ppbv was observed at DRI in 2013, when using a Membrana 160 Micromodule degasser for analysing the Tunu13, D4, and a fraction of the NEEM-SC cores. Lower blank values, ranging 4.1±1.2 to 12.6±4.4 ppbv were observed when using an Idex Transfer-Line degasser (cores PLACE, NGRIP, a fraction of NEEM-SC).

Internal precision and stability
Internal precision and stability of gas-CFA measurements can be evaluated by Allan variance tests (Allan, 1966) applied to the 165 calibration loop dataset. Laser spectrometers produce many measurements per minute, with a data acquisition rate of 6 Hz for our OF-CEAS instrument. Observed optimal integration time (i.e., time of lowest Allan deviation) was determined for each analytical campaign, and is now larger than 500 s for DRI, and 1000 s for IGE CFA setups (SI), which is much longer than the 5 s reported by Faïn et al. (2014) and testament to increased stability. However, to maximize the depth resolution, CO data can be averaged over shorter integration time (IT). There is a trade-off between internal precision and IT. Abrupt, non-climatic 170 signals in gas records, such as in situ CO production (Faïn et al., 2014), cannot be fully resolved without reducing the IT.
Internal precision, defined as twice the Allan deviation at chosen integration times, ranged from 1.2 to 1.6 ppbv depending on the analytical campaign (SI :::: Sect. ::: 1.6 : and Table S2).

External precision
External precision of the continuous CO measurements (i.e., including all sources of errors or bias) was investigated by melting 175 replicate ice sticks on different days during each analytical campaign. Specifically, we defined the external precision as the pooled standard deviation calculated on the differences of CO concentrations from main and replicate analysed ice sticks, averaging continuous CO data over few cm long intervals :: (SI ::::: Sect. :::: 1.7). This approach estimates an external precision of 5.7 ppbv for the DRI gas-CFA setup in 2011 (NEEM-2011-S1 campaign, 18 m long replicated section, using 9 cm intervals). We reproduced this analysis for the 2015 DRI campaign, and established an external precision of 6.6 ppbv (19.6 m long replicated 180 section, 10-cm long intervals, Fig. S3). With a similar methodology, but a shorter replicate section analysed, we evaluated the external precision of PLACE continuous CO measurements conducted at IGE in 2017 to 3.3 ppbv (2.3 m long replicated section, 1-cm intervals, Fig. S3). No replicate ice core sections were available and measured during the 2013 DRI campaign.
For the 2013 campaign, we consider the external precision to be 5.7 ppbv for the ice cores analyzed in 2013 with a Micromodule and 6.6 ppbv using the Idex Transfer-Line degasser, by extrapolating results from similar CFA setups. 185

Absolute calibration and accuracy
CFA-based gas records must be corrected for under-recovery of gases dissolved in the water stream to obtain absolute values on the WMO-CO X2014A calibration scale. The magnitude of gas dissolution is monitored using the calibration loop (Sect. 2.2.2). In this study, a 73 m long ice core (PLACE) was analysed on two different CFA setups (at IGE and DRI), for CO but also for CH 4 mixing ratios. This gas CFA laboratory intercomparison revealed that the fraction of CH 4 not recovered at IGE 190 was larger than at DRI (14% and 10%, respectively). Note that the gas extraction unit and OF-CEAS detection instrument were identical, and only the melter and lines upstream of the Idex Transfer-line degasser had different geometries (Sect. 2.2.3).
This intercomparison demonstrates that the solubility-based calibrations are required to transfer gas CFA datasets into absolute concentration scales, and are dependent on the geometry and operation of a CFA setup. This observation was further confirmed at IGE by measuring the CH 4 mixing ratio at IGE when melting replicate ice sticks with varying CFA system geometries (eg., 195 length of lines) and operation of the CFA setup (e.g., melting speed, pumping rate). For each analytical campaign, the geometry and operational protocol of the CFA setup were thus kept unchanged during the measurements. However, these observations rule out universal modelling of the preferential CO dissolution by using Henry's Law coefficient, and demonstrate that CFA setups are dynamic systems, not at solubility equilibrium.
To further evaluate the robustness of the CFA procedure and this absolute calibration approach, five discrete samples of the

Signal Smoothing
The mixing of gases and meltwater during the sample transfer from the melt head to the laser spectrometer induces a CFA experimental smoothing of the signal. The extent of the CFA-based damping was determined for each analytical campaign by performing step tests, i.e., switches between two synthetic mixtures of degassed DI water and synthetic air standards of 220 different CH 4 concentrations (Stowasser et al., 2012, see SI :::: Sect. ::: 1.9). A cut-off wavelength can be defined as the wavelength of a sine signal experiencing a 50% attenuation in amplitude. Cut-off wavelengths ::::: ranged : from 15.0 cm during the DRI 2013 campaign (when using Membrana Micromodule degasser) to 1.6 cm during IGE 2017 campaign (Table S2). Lower cut-off wavelengths were obtained by limiting system dead volume, most notably by reducing the debubbler internal volume and by using an optimized Idex Transfer-Line degasser instead of Membrana Micromodule degasser (1 and 5 ml internal volume, 225 respectively)

Data processing
CO and ice chemistry data collected at DRI were mapped onto depth scales using high-resolution (0.1-0.5 Hz acquisition rate) liquid conductivity data and time-depth relationships recorded by system operators. A constant melt rate for each metre length of core is assumed. DRI depth scale uncertainties are estimated to be ±2 cm (2σ : 1 ::: cm :::: (1σ).

230
Occasional entry of ambient air into the analytical system as breaks in the core were encountered caused contamination. The OF-CEAS spectrometer simultaneously measures CO and CH 4 mixing ratios, and such contamination was characterised by a sharp increase in CH 4 concentration (ca. 1900 ppbv) followed by an exponential decrease. Data were manually screened for these ambient air contamination events.
Uncertainty on CO mixing ratios is established as 1σ and calculated for each calibrated CO dataset. It combined uncertainties 235 evaluated specifically for each analytical setup on CO blanks, solubility calibration factors, and internal precision of CO CFA measurements (see SI and previous sections).
Finally, we used high resolution water isotope datasets to match DRI and IGE PLACE depth scales. Water isotopic ratios were measured simultaneously with gas concentrations in both laboratories using laser spectroscopy (Gkinis et al., 2011;Maselli et al., 2013). To investigate if elevated or highly variable carbon monoxide levels observed in Greenland ice core could originate from the melting process (i.e., "in extractu" CO production from trace organics in the ice), CO concentrations were determined in tests with 20 discrete ice and firn samples (sized between 272 -1089 g, average size of 690 g). The experimental protocol is reported in detail in the SI. Briefly, ice samples collected below close-off depth (fully closed porosity) are completely grated to 245 fine powder to remove trapped air. Firn samples do not contain a significant amount of trapped air; however some firn samples were still grated to powder to verify that the grating process did not increase the [CO] blank. The samples are then introduced in a glass vessel, evacuated and flushed with ultrapure air overnight. Gas standard of known concentration is then introduced within the vessel. Consequently, for both firn and ice samples CO mixing ratio within the vessel prior to melting is equal to that of the gas standard introduced in the vessel. If CO in extractu production occurs during melting, it should cause an increase in 250 CO mixing ratio above this known value.
The CO concentration in the gas available in the vessel headspace is determined by gas chromatography combined with a mercuric oxide reduction detector (RGD; Peak Performer 1 from Peak Laboratories), at three different stages: before melting, during melting, and after melting. CO concentrations are calibrated using four synthetic air gas standards with nominal CO concentrations ranging from 50 to 500 ppbv, and data are reported on the WMO-CO X2014A scale.

255
Three large diameter (27 cm Antarctica ice, and (v) gas-free ice made from degassed Milli-Q 18.2MΩ.

Chemistry data
During the two DRI analytical campaigns, melted ice core samples were analysed continuously by inductively coupled plasma 265 mass spectrometry (ICP-MS) and CFA for chemical species. These analytical methods have been reported previously (Mc-Connell and Edwards, 2008;McConnell et al., 2007McConnell et al., , 2002. High-resolution measurements of total organic carbon (TOC) were obtained at DRI by coupling a Sievers 900 TOC analyzer to the DRI ice core melter (Legrand et al., 2016).

High frequency non atmospheric CO signals
The true extent of high CO variability may be masked, however, by CFA analytical smoothing. Figure 2 (upper panel) reports a 1.5 m long section of the PLACE ice-core, analysed with both DRI and IGE CFA setups. The higher CFA setup resolution 300 at IGE (Sect. 2.2.7) reveals more features, but may still miss even higher frequency variability and dampen extrema values.
Similarly, a CO record from a core drilled at a relatively low snow accumulation site will appear more smoothed out by CFA analysis than one from a relatively high snow accumulation rate site. Therefore, in this study we do not directly compare the amplitude of CO variability between D4 and Tunu13, which exhibit a factor of four difference in snow accumulation rate.
High frequency variability within CO records can be quantitatively characterized using MAD (Median Average ::::::: Absolute : De-305 viation). MAD was calculated for all CO records collected at DRI every 4 yrs over a moving window of 15 yrs (similarly to baseline signals, see Sect. 3.1.1) :: are :::::::: reported :: on :::: Fig. :::: S18. We observe for the 1700-1950 CE time period that all :::::: NEEM, clearly indicating that high variability in CO is more pronounced for older ice. These increases are significant, with MAD values reaching ∼ 60 ppbv for the deeper sections of the analyzed NEEM and NGRIP ice-core sections.
3.2 Do drilling or analytical processes cause CO production?

315
The observation of high frequency CO variability along the NEEM-2011-S1 core led Faïn et al. (2014) to discuss if such patterns could be driven by the processes involved in sample collection and analysis. Here we consider our new data to assess in more detail if such processes can drive CO production.

Ice core drilling
We observe that abrupt CO spikes of similar magnitude occur in ice-cores extracted by dry drillings (e.g., PLACE, D4, Tunu13) 320 and when a drilling fluid was used (e.g., NEEM, NGRIP) ( Fig. 1). This clearly rules out contamination from the drilling fluid as a cause for high frequency CO variability, confirming the previous conclusion from Faïn et al. (2014). Haan et al. (2001) observed CO production within an alpine snowpack in daylight. Drillings are usually conducted at highaltitude sites in summer where UV radiation levels are high. One could thus hypothesise that CO is photochemically produced within ice-cores just after drilling, when ice-cores are handled at the drilling site before being packed in ice-core boxes. During 325 the PLACE and NEEM-SC drillings, specific care was taken to never expose directly freshly drilled cores to sunlight. For PLACE, a shaded area was set up for core processing. NEEM 2011-S1 and SC core drillings were conducted in a tent and a trench, respectively. In spite of that, all these ice-cores reveal CO spikes. In contrast, no specific care was taken during the historical Eurocore drilling when handling cores, with respect to sunlight exposure, and Haan and Raynaud (1998) report that this archive does not exhibit highly variable CO concentration in the upper section. Based on this, it appears that brief direct 330 exposure to sunlight after drilling is unlikely to be the cause of the large CO spikes 3.2.2 Impact of long-term core storage By comparing the storage time of different ice-cores (storage without light exposure) prior to CO analysis with the same CFA setup (NEEM-2011-S1 and D4), Faïn et al. (2014) concluded that CO production in ice during core storage was unlikely. We extend this observation by comparing central Greenland archives: NGRIP and PLACE records reveal similar abrupt CO spikes 335 and CO baselines ( Fig. 1) while PLACE and NGRIP cores were analysed 3 months and 15 years, respectively, after drilling.
Furthermore, we measured at IGE a 10 m long replicate section of the PLACE core nearly two years after the main analytical campaign (respectively December 2018 and February 2017) and no change in the CO values was observed (Fig. S12 ::: S14).

Impact of CFA process
In 2013 at DRI, a 4 m long replicate section of the D4 ice-core was melted with all laboratories kept in darkness. Similar in situ 340 CO production patterns were observed under light and dark melting conditions. These results, already reported by Faïn et al.
(2014, Fig. S2), reveal that CO variability observed along the D4 core are not related to photochemically driven CO production occurring in the sampling lines between the melter head and the CFA extraction box. To further evaluate if the CFA analytical system could somehow induce CO production within the melted ice, a melted D4 sample was collected downstream of the degassing membrane and re-circulated in the calibration loop mode (Sect. 2.2.2) for one hour as degassed blank water. We 345 observed similar CO levels to the deionised water blank.

3.2.4
Results of investigation of possible rapid CO production from trace organics in the ice during melting As described above, CO analyses of 20 discrete ice and firn core samples (including the PLACE core) were carried out at the University of Rochester, with the goals of investigating the possibility of in extractu production of CO from trace organics in the ice during melting of ice core samples (Sect. 2.3 and SI :::: Sect. :::: 1.11, Table S4).

350
All of the different ice and firn sample types, including samples of different mass and a gas-free ice sample, showed a consistent excess CO growth occurring during the melting step, with excess CO increasing by 6.4 ppbv on average (Fig. 3) -a much smaller CO enhancement than the spikes observed in the CFA records. Since the excess CO produced appears stable after melting, and consistent between the sample types, the excess CO observed is likely a systematic extraction system blank.
If the CO production mechanism involves organic compounds present in the ice lattice, it would be expected that variations in T i m e a f t e r m e l t s t a r t s ( m i n ) G a s F r e e I c e C e n t r a l G r e e n l a n d F i r n C e n t r a l G r e e n l a n d   Table S4. Error bars represent one standard deviation for all sample runs of a given type and points without error bars represent single measurements . the concentrations of organics from samples collected at different sites / from different depth levels, as well as differences in the melted ice mass would result in significant differences in excess CO. Some of the PLACE ice samples were specifically selected for having exhibited low and stable [CO] on the continuous CFA record while others were selected for exhibiting elevated and spiky [CO] behavior suggesting the possibility for instantaneous CO production during melting (Fig. 3, Table S4).
The other samples, the two Taylor Glacier and gas-free ice were also selected as they should have very different trace organic 360 loadings to the Greenland ice, especially the gas-free ice sample which have no significant organic loading (< 2 ppbC). Further, replicate firn samples of different mass would be expected to have a total organic loading that scales with the firn mass. Because there are no significant differences observed in the excess CO produced from these very different sample types (Fig. 3), these results indicate that the melt-extraction process itself does not result in significant CO production from trace organics found within the ice samples. Instead, these results lend further support to the idea that the elevated and highly variable [CO] values 365 observed in Greenland ice samples are due to excess CO produced from in situ production within the ice itself (Faïn et al., 2014;Haan and Raynaud, 1998).

Atmospheric CO history retrieved from Greenland ice cores
We now explore if past atmospheric concentrations of CO can be extracted from the low frequency variability of the CO icerecords' baselines ( Fig. 1). The key question is: does the in situ production present in all five cores still impact the 5 th percentile 370 of data that we adopt as a baseline?
These are interesting cores to compare because they have very different snow accumulation rates: Tunu13 has the lowest accumulation rate, which also varies over time between 8 to 13 cm weq yr −1 , while PLACE is more typical of a Late Holocene

380
Greenland core with a stable accumulation rate of ∼20 cm weq yr −1 .
Typical CO, TOC, and NH + 4 variations along the Tunu13 core appear different to those at PLACE (Fig. 4). The seasonal 390 cycle in NH + 4 is still detectable but the CO and TOC :::: TOC :::: and ::: CO records show much less co-variation than at PLACE.

TOC patterns and in situ CO production [REMOVED SECTION]
The high resolution TOC measurements collected along the PLACE core exhibit a stable TOC background value of 13.5 ppbC (baseline defined as the 5 th percentile), with 17 spikes above 40 ppbC (70 m long record, Fig. S18). These spikes are 455 often co-located in depth with increases in CO concentration (e.g., Fig. 2) as observed previously on the NEEM-2011-S1 core (Faïn et al., 2014), suggesting that in situ CO production may be related to the organic carbon (OC) availability. TOC levels observed at Tunu13 are lower, with a baseline value of 7.9 ppbC (1σ) and only 10 spikes above 40 ppbC for a longer record (143 m long record, Fig. S19). However, mean TOC (10 meters window averaging) exhibits a low frequency variability with values ranging from 4 to 12 ppbC (Fig. 5). This low frequency change is correlated to change of the snow accumulation rate and snow accumulation rate remains unchanged, although it becomes steeper below 160 m depth (Fig. S20). Such correlation is not observed at PLACE where the snow accumulation is constant, with annual mean values ranging 20-21.5 cm weq yr −1 . Our results indicate that the fate of TOC deposited at Tunu13 is related to the rate of snow accumulation at the surface. Both Tunu13 and PLACE datasets suggest that the lower the accumulation rate, the lower the TOC baseline concentrations in the deep ice.

465
A snowpit study of recent snow layers deposited at Summit (Hagler et al., 2007) showed that a large fraction of water soluble and insoluble organic carbon may be lost at the surface during post-deposition processes, such as photochemical reactions.
We hypothesize that such processes may be enhanced at lower accumulation rates when the snow layers get exposed to UV radiation for longer. However, such processes may not be the only ones at play, as suggested by investigating the relationships between TOC and ammonium. Forest fire debris reaching Central Greenland snow mainly consists of ammonium formate, 470 and organic carbon is mostly made of formate at Summit (Legrand et al., 2016). Legrand and De Angelis (1996) reported a formate-to-ammonium molar ratio close to unity in elevated ammonium events recorded at Summit both during the last 200 years and the Holocene. The relationship between TOC and ammonium at PLACE (Fig. S21) exhibits a similar slope close to unity. A completely different pattern is observed at Tunu13 (Fig. S22), with unexpectedly low TOC observed along ammonium peaks. We suggest that ammonium formate has been remobilized after deposition in the snow at sites with snow accumulation 475 rates lower than 12 cm weq yr −1 . While ammonium remains when deposited due to acidic snow layers, formate evolves to a different form, presumably gaseous formic acid that, in contrast to ammonium, can escape from neighbouring acidic snow layers and possibly partly be released into the free atmosphere above the snowpack. Such a process could contribute to explain lower TOC levels observed under lower snow accumulation rates, and also the smoothed shape of the TOC peaks at Tunu13 (Fig. 4).

TOC patterns and in situ CO production [REMOVED SECTION]
In situ CO production within Greenland ice archives is very often co-located in ice layers where TOC levels are high. The importance of organic carbon in the processes driving CO in situ production is supported by the long-term averaged mean CO record (running windows of 10 m) which is significantly correlated to mean TOC at Tunu13 (r 2 = 0.56, p<0.01, Fig. 5 (20 and 7 ppbC, respectively). This may suggest that the required amount of OC to lead to significant in situ CO production is already reached at the Tunu13 site. While the mean CO mixing ratio is significantly correlated to the mean TOC concentrations 490 in ice at Tunu13, this is not the case for baseline CO level which does not exhibit significant correlations with mean or baseline TOC concentrations. However, similarities in trends and patterns can be seen when comparing Tunu13 baseline CO and TOC ( Fig. ??). The larger analytical smoothing impacting the Tunu13 CO record means that some of the CO baseline signal likely incorporates in situ produced CO from spring/summer ice layers. We can not rule out, however, that a redistribution of organic carbon along depth driven by OC post-deposition process (shifting the ammonium-formate equilibrium with the gas phase) Running averages (Continuous (high-resolution) and baseline (5 th percentile calculated over a 1 m window) CO (red) and TOC (blue) records for the Tunu13 ice core.

Potential of low snow accumulation sites for long-term CO reconstruction [REMOVED SECTION]
There are signs that over the long term (i.e., 500 yrs and more), low accumulation sites such as Tunu13 are promising for 500 reconstruction of paleoatmospheric CO records. For example, the deepest section of the Tunu13 CO record does not show an increase in MAD (Sect. 3.1) and exhibits a relatively stable baseline. This is different from higher accumulation sites such as NGRIP or NEEM, where MAD increases with depth, and a positive trend in CObaseline with depth is clearly observable (Fig.   S15). The chemical processes involved in CO in situ production are likely multiple and complex and the oxidants involved remain unknown. Nevertheless, we observe that H 2 O 2 , a potential oxidant, exhibits lower levels when accumulation is lower 505 (i.e., Tunu13 vs PLACE, Fig. S23).

Comparison of PLACE and
To further investigate this apparent offset between our new Greenland CO records and the Eurocore CO record of Haan et al., we analysed six 50 cm-long sections of the 1989 Eurocore ice core with the gas CFA IGE setup in January 2019, at depths 535 spanning 84 to 153 m. A fraction of the Eurocore core has been stored as an archive for 30 yrs at -20°C. We were able to sample the Eurocore archive at exactly the same depths that Haan et al. did in the 1990s (the gaps left by the discrete sampling conducted by Haan and colleagues could be seen). Figure 6 shows Eurocore data (both from Haan et al. and our new CFA dataset), along with the PLACE IGE CFA record: to do so, an offset of +6 m was applied to the Eurocore depth scale, so as to take into account 26 yrs of snow accumulation that occurred between the Eurocore and PLACE drillings. Our CFA-540 based CO analyses reveal high variability along the Eurocore archive. We systematically observed higher CO mixing ratios than previously reported by (Haan and Raynaud, 1998) (Fig. 6), while they show an excellent agreement with the CFA CO dataset (Fig. S7).
To summarize, we have not been able to reconcile the historical Eurocore dataset of Haan et al. with a new Greenland ice-core CO dataset (PLACE). This is despite analysis of remaining Eurocore ice from depths co-located with the published Haan data and the fact that PLACE ice core was collected just 1 km away from Eurocore. The reasons for the relatively low and stable 550 CO measurements reported in the 1990s for the Eurocore ice core remain unknown and our results bring the validity of the Haan et al. (1996,1995) results into question.
From 1875 to 1950 CE, the records indicate a monotonic increase from 100-120 ppbv to 135-150 ppbv (i.e., a mean rate of increase of ∼ 0.3 ppbv yr −1 ). Overall, Fig. 1 reports a ∼ 30% increase in CO concentration at high latitudes of the northern hemisphere atmosphere from pre-industrial to 1950 CE. Interestingly, in 1835 CE, PLACE, NGRIP, NEEM, and D4 baselines 565 exhibit a common maximum in CO concentration but at a level just slightly higher than concentrations observed during the 1700-1820 CE period.
A multisite composite CO baseline spanning 1700-1957 CE was generated by considering D4, NEEM, NGRIP, and PLACE records and excluding Tunu13. For consistency, we only consider CO records collected on the same CFA gas setup (DRI).
The multisite composite, reported in Fig. 7, was obtained by averaging CO baselines values interpolated on a common gas age  The fact that the CO baseline records from the four different Greenland ice cores are so consistent gives us confidence that our multisite reconstruction provides a reliable reconstruction of past atmospheric changes. An additional argument for that conclusion comes from the good overlap between our multisite reconstruction and CO data from firn air measurements con-575 ducted at Summit, Greenland (see section 3.5.2). However, we cannot exclude the possibility that past atmospheric CO levels could have been lower than shown in Fig. 7 and recommend that our multisite average should be taken as an upper bound of past atmospheric CO burden. If in situ CO production drives the high frequency variability seen throughout the CFA records of these four ice cores, as we hypothesize (see Sect. 3.2), it is impossible to exclude the possibility that minimum CO levels may also be slightly affected (Sect 3.3).
Data availability. The high resolution carbon monoxide datasets will be made available on the World Data Center for Paleoclimatology.
Author contributions. This scientific project was designed by XF, JC, VVP, RR, EB, and TB. The high resolution carbon monoxide mea-    Table S1. Information on ice and gas age scales for the ice cores featured in this study 9

Differences and similarities between DRI and IGE CFA setups
1 Specific details of analytical setups are reported in Table S2. Small operational differences exist 2 between the laboratory setups at DRI and IGE. First, DRI and IGE melter geometries are different, 3 and described by McConnell et al. (2002) for DRI, and Bigler et al. (2011) for IGE. Second, the 4 melting rate was higher at DRI compared to IGE ( [5][6] cm min -1 and 4 cm min -1 , respectively). 5 Finally, DRI and IGE debubblers have different geometry. DRI debubbler (Plexiglas-made) allows 6 for ultra-high resolution of liquid phase analyses, but results in a larger mixing of the gas phase 7 sample compared to the IGE debubbler (glass-made) specifically designed for gas measurement. 8 The same optical feedback cavity enhanced absorption spectrometer was used to analyse carbon 9 monoxide and methane at both DRI and IGE and during the all analytical campaigns (2013,2015,10 and 2017). The same instrument was used by Faïn et al. (2014) for measurements of the NEEM-11 2011-S1 core in 2011 at DRI. 12 Two different gas extraction units were used (i) in 2013 (DRI), and (ii) both in 2015 at DRI and 13 2017 at IGE, respectively. A gas extraction unit includes a degasser, a pressure transducer, a 6-14 ways valve, and a homemade Nafion dryer (Fig.1 of Rhodes et al., 2013). Improvements on the 15 design and operation of this unit largely explain the improved stability of the [CO] CFA signal, as 16 shown by the much larger optimal IT values observed for the 2015 and 2017 campaigns (Table  17   S2 Table S2. CFA setup specifications and performances of the different analytical campaigns. 1 The optical response of OF-CEAS (based on absorption) to gaseous concentrations of samples 2 circulating in the cavity of the spectrometer is highly linear . Such linearity has been 3 previously reported for concentrations ranging over more than three decades for methane (e.g., 4 Fourteau et al., 2017). 5

Calibration of SARA analyser
Similarly, we observed a highly linear relationship of our OF-CEAS spectrometer when measuring 6 CO mixing ratios over a [0-100 ppbv] range. Such linearity was demonstrated by a dilution 7 experiment conducted by mixing a 100 ppbv standard gases and CO -Air Zero with two MFC 8 MKS recently calibrated (concentrations ranging [0-100] ppbv could thus be generated). For each analytical campaign, the OF-CEAS spectrometer was carefully calibrated on dry gas by 1 direct injection of a synthetic standard gas (Scott Marin, artificial gas mixtures) precisely calibrated 2 onto the NOAA/WMO X2014 scale, with a concentration close or above 100 ppbv (Table S3). 3 Knowing the excellent linearity of the instrument, a calibration factor was established as the ratio 4 of the NOAA-certified and the measured CO concentrations for this single standard gas. This is 5 consistent with the fact that the zero of the spectral measurements is intrinsically accurate. 6 7  how to quantify exactly such procedural blanks remains complex. Two approaches investigate 15 gas-CFA procedural blanks: single-standard measurements, or multi-standards calibration. 16 To evaluate CO blank with single-standard measurements approach, deionized (DI) water is 17 bubbled with a specific gas standard for at least 12 hours in a 2 liters reservoir, and this water is 18 later mixed with the same gas standard in the calibration loop (CL, section 2.2.2 of manuscript). 19 The positive offset between direct measurements of the dry gas standard with the spectrometer 20 and the CL measurements as described previously will directly indicate the procedural blank. 21 Single-standard measurements were conducted to evaluate CO blank of the IGE CFA setup. An 22 example of such experiment is as follow: DI water bubbled during 18 hours with a 58.7±0.3 ppbv 23 (NOAA/WMO X2014 scale) CO standard gas. The same standard gas was then mixed with this 24

25
DI water within the IGE calibration Loop, and circulated toward the degasser. CO mixing ratios of 1 62.8±0.5 ppbv is measured. A procedural CO blank of 4.1±0.6 ppbv is thus observed. Multi-standard calibration loop datasets have been widely used before for methane (e.g., Rhodes 8 et al., 2016;Fourteau et al., 2017) and are obtained by mixing a range of standard gases of known 9 concentration with He-degassed DI water. The intercept observed when plotting expected versus 1 CL-determined concentrations reflects both the procedural blank and the fraction of gas remaining 2 in DI water after He-bubbling. Figure S4 reports multi-standards CO calibration obtained at DRI 3 in 2013 with both Micromodule and Transfer Line degassers, and the intercepts of 35 ppbv and 4 12.6 ppbv are the CO blank levels observed for these two specific CFA configurations. We applied 5 the multi-standards calibrations loop approach to evaluate CO blank for all continuous CO 6 datasets collected in this study (Table S2). At IGE, the multi-standards calibration provided an 7 Intercept for [CO] of 4.1 ppbv, in good agreement with the single-standard approach described 8 before. This good agreement demonstrates that the multi-standards calibration is a reliable 9 approach to evaluate CFA blank for gas of low solubility such as CO, i.e. the fraction of CO 10 remaining in DI water after He-bubbling is negligible. Estimation of the uncertainty of CO blanks 11 are reported on Table S2 and are based on repeated calibration loop measurements throughout 12 the campaigns (1σ). Also, any material (e.g., line, connections) should be carefully tested to assess potential CO 3 contamination. We have found that TEFZL is appropriate for CO analyses. Finally, CO blanks are 4 kept low by continuously running DI water through system lines during measurements periods. 5 6 1.6. Internal precision of CO CFA analyses 7 Internal precision and stability of gas-CFA measurements were evaluated by Allan variance tests 8 (Allan, 1966)

External precision of CO CFA analyses
6 External precision of the continuous CO measurements (i.e., including all sources of errors or 7 bias) can be investigated by melting replicate ice sticks on different days on a gas-CFA setup. 8 Specifically, we define the external precision as the pooled standard deviations calculated on the 9 differences of CO concentrations from main and replicate analysed ice sticks, averaging 10 continuous CO data over few cm long intervals. Specifically, after binning the main (M) data and 11 replicate (R) data into these few cm long intervals, we obtain n duplicate measurements. The 12 pooled standard deviation is then calculated as √(( For the DRI gas-CFA setup, an external precision of 5.7 ppbv can be extracted from a 18 m long 14 replicate section of the NEEM-2011-S1 core analyzed in 2011 (Faïn et al., 2014) Fig. S3. CO concentrations measured on both main and replicate ice sticks were averaged by 8 binning over 196 10-cm long (resp., over 227 1-cm long) intervals at DRI (resp., at IGE). Figure  9 S6 reports an excellent correlation (r2 = 0.99; p<0.01) between averaged CO from main cuts 10 versus averaged CO from replicate cuts for both CFA setup. Finally, external precision of 6.6 ppbv 11 (resp., 3.3 ppbv) was calculated for the 2015 DRI campaign (resp., 2017 IGE campaign). In this 12 study, we apply the external precision value established for the 2015 DRI campaign is applied to 13 dataset collected at DRI in 2013 with a Transfer-Line degasser (i.e., a fraction of NEEM-SC and 14 NGRIP dataset). 15 1 Figure S6. CO mixing ratios measured on both main and replicate PLACE ice sticks at DRI (a) and IGE 2 (b), after averaging data over 10-cm long (resp., over 1-cm long) intervals at DRI (resp., at IGE). 3 4 We are aware that longer replicate sections, spanning depths over the entire record, such as 5 conducted at DRI in 2015 (Fig. S3), can provide a more solid dataset for evaluating external 6 precision. We note that only the DRI PLACE dataset was used to extract atmospheric information 7 (section 3 of manuscript), notably because a more robust evaluation of the external precision was 1 available. However, the IGE dataset demonstrates that excellent precision of continuous CO 2 measurements is not limited to one specific setup, but rather achievable for all CFA laboratories. 3 4 1.8. Internal CO CFA calibration: improving accuracy 5 A fraction, or all, of gases dissolved in the water CFA sample flow are not recovered by the 6 degassers used in this study. CO solubility is higher than N2 or O2 solubility. This preferential 7 dissolution of CO in comparison to N2 and O2 results in underestimated CO mixing ratios when 8 detected at the gas outlet of the degasser. This underestimation of CO mixing ratios need to be 9 accurately quantified so as to provide an absolute calibration of CO continuous dataset. Henry's 10 Law theory suggests that temperature and gas/liquid ratio of the melted sample are two important 11 parameters influencing gas dissolution within the sample flow. Gas/liquid ratio of the melted 12 sample is likely related to the air content of the ice analysed, although this relationship may not 13 be straightforward (e.g., one cannot exclude losses of gas at the melter), and may differ with 14 melter geometry. 15

Rationale for an internal CFA calibration 16
In this study, we use the calibration loop so as to calibrate internally (i.e., with CFA data) 17 continuous CO datasets. This choice relies on the following observations: 18 • The amount of gas dissolved in the water stream and not recovered by the degasser (Idex 19 Transfer-Line or Membrana Micromodule) is dependent on the CFA setup, i.e. dependent 20 on its geometry, components, operation. Consequently, a unique theoretical framework 21 cannot be applied to all dataset. Indeed, in this study a 73 m long ice core (PLACE) was 22 analysed on two different CFA setups (at IGE and DRI), for CO but also for methane 23 mixing ratios. This gas CFA laboratory intercomparison revealed that the fraction of 24 methane not recovered at IGE was larger than a DRI (10% and 13%, resp.). The gas 25 extraction unit (including an Idex Transfer-line degasser) and OFCEAS detection were 26 identical at DRI and IGE, pointing out that the melter and lines upstream the degasser can 27 result in different gas dissolution ratios. 28 • To test further the influence of the CFA setup geometry and operation on gas dissolution 29 patterns, a series of 16 replicate 1m long ice sticks, all cut from the same Blue Ice Drill Specifically, the influence of following parameters on the preferential methane dissolution 3 were tested: melting speed (varying from 2.8 to 6 cm min -1 ), sample flow (ranging 12 to 22 4 ml min -1 ), melter geometry (geometry reported by Bigler et al., 2011, with inner sections 5 of 26 x 26 mm, or 24 x 24 mm), line lengths (adding 4 m long line sections in the system). 6 These tests are specific to the IGE setup, but they show: (i) no influence of the melting 7 speed, neither of the melter geometry, neither of the sample flow, on the preferential 8 methane dissolution. However, they reveal varying methane losses when changing line 9 geometry (up to 5% in mixing ratio). 10 Overall, these results show that CFA setups are dynamic systems, not at solubility equilibrium. 11

Principle of calibration 12
We hypothesize that CO and methane dissolution follow the same physical laws: consequently, 13 if a calibration loop (CL) is able to reproduce methane preferential dissolution, it should also 14 reproduce CO losses related to dissolution. Calibration loop is described in sections 2.2.2 of 15 manuscript and 1.5 of SI; briefly, a 10:90 mixture of synthetic air and degassed deionized (DI) 16 water can be introduced into the system via a 4-port valve located directly after the melterhead. 17 The water is sourced from a 2 L reservoir degassed by constantly bubbling ultra-pure helium 18 through it, and the fraction of He dissolved in DI water is later neglected. The air-water mixture 19 follows the same path through the system as the ice core sample before being analysed by the 20 laser spectrometer. We verified that no CO or CH4 specific fractionations occur at the degasser, 21 which would invalidate the previous hypothesis. To do so, we run consecutively the calibration 22 loop in two configurations: (i) with a degasser, and (ii) with an open split debubbler where gases 23 are transferred to the OF-CEAS analyzer without going through a membrane. Same CH4 and CO 24 mixing ratios were observed in both configurations. Finally, considering that gas dissolution within 25 the CFA setup is likely to vary with ice air content (AC), we tried to determine a specific calibration 26 for each ice core analyzed. 27 We first needed an estimation of methane losses related to dissolution, independent of the 28 calibration loop, from each ice core analyzed. For the archives Tunu13, D4 and NEEM, 29 preferential methane dissolution was already known (Rhodes et al., 2013(Rhodes et al., , 2016. For the PLACE 30 and NGRIP archives, it was evaluated by direct comparison of CFA dataset with GISP2 discrete 31 record . were analysed discreetly, as this analytical capability was designed only in late 2019, and applied 1 in priority to investigate if high frequency CO spikes observed in Greenland continuous CO 2 records could be driven by in-extractu, production as described in Sect.  production occurs during the melting analytical process). However, the depths of the five discrete 4 samples described here ranged the entire core, from 84 to 146 m depth. 5 A direct comparison of discrete and continuous dataset is reported on Fig. S7. The standard 6 deviation of the differences between discrete and CFA-based CO measurements is 1.8 ppbv and 7 2.8 ppbv for the IGE (n=5) and DRI (n=4), respectively. The average offset between discrete and 8 CFA-based CO measurements is 1.3 ppbv and 2.1 ppbv for the IGE (n=5) and DRI (n=4), 9 respectively. This excellent agreement supports our approach to internally calibrate continuous 10 CO, and demonstrates that continuous [CO] dataset can now be accurate, which is an important 11 improvement since the NEEM-2011-S1 measurements (Faïn et al., 2014). 12 13 Figure S7. Comparison  Due to mixing and dead volumes, the CFA system introduces a smoothing of the gas signals. 4 Using the method described in Stowasser et al. (2012), we measured during each analytical 5 campaign the switch between two mixes of deionized water and standard gases ( Fig. S8-11, left 6 panels). It allowed us to determine step responses of the CFA systems, which is not instantaneous 7 but spreads over time. During the 2013 DRI campaign, both Idex Transfer-Line and Membrana 8 Micromodule degassers were used (Table S2), and step responses were determined for each 9 degasser configuration. 10  Figure S10. Same as Figure S8 for the DRI CFA system operated with a Transfer Line (Idex) degasser in 6 2015 for analyzing the PLACE ice core. 7 8 Figure S11. Same as Figure S8 for the IGE CFA system operated with a Transfer Line (Idex) degasser in 2 2017 for analyzing the PLACE ice core. 3 4 Each observed step response was fitted using the cumulative density function of a log-normal 5 distribution. These log-normal distributions ( Fig. S8-11, middle panels) can be considered as the 6 Green's functions, or impulse responses of the CFA systems. Finally, the Green's functions can 7 be used to derive the frequency responses of the systems using Fourier transforms, that is to say 8 the attenuation factors (also referred to as gains) experienced by a sine signal depending on its 9 frequency/period/wavelength. A cut-off wavelength can be defined as the wavelength of a sine 10 signal experiencing a 50% attenuation in amplitude. It is important to note that this cut-off is 11 defined for sine signals, and therefore cannot be directly applied to other types of signals. For 12 instance, Fourteau et al. (2017) reports that a square signal will be less attenuated than a sine 13 signal. 14 During the 2013 DRI campaign, we observed with a melting rate of 5.5 cm min −1 cut-off 15 wavelengths of 15.0 cm and 14.2 cm, when using a Micromodule or a Transfer-Line degasser, 16 respectively. A Micromodule was used for the D4, Tunu13, part of NEEM-SC, cores, and a 17 Transfer-Line degasser was used for NGRIP, and a fraction of NEEM-SC (Table S2). During the 18 2015 DRI campaign, we used an Idex Transfer-Line degasser optimized, i.e. shortened to reduce 19 its internal volume to only 1 mL. Consequently, we observed with a melting rate of 5.5 cm min −1 20 a cut-off wavelengths of 9.3 cm. During the 2017 IGE campaign, we used the same shortened 21 Transfer-Line degasser, but the entire gas setup was optimized specifically for gas measurements 22 (including the use of a low dead volume glass debubbler): a cut-off wavelength of 1.6 cm was 23 observed. Similar improvements were investigated on the DRI CFA setup in 2015, but could not 1 be fully implemented because they would have decreased the resolution of the liquid phase 2 analyses conducted simultaneously to CO continuous measurements. The Tunu13 ice core was analysed with the 2013 DRI CFA setup (Micromodule degasser), which 7 exhibit the largest response time in the gas phase (Table S2, Fig. S8). This effect is compounded 8 by the low snow accumulation rate of Tunu13.
Step tests (Sect. SI 1.9.1) suggest that for CO a 9 60% attenuation would be expected for a sine signal of 10 cm wavelength, a depth interval 10 representing >1 yr. 11 To further investigate how analytical smoothing can alter the CO 5 th baseline levels of the Tunu13 12 archive, the following modeling was conducted: 13 -A synthetic signal representing the original Tunu13 CO signal is produced as follow: the 14 baseline is set to a constant 80 ppbv value, with abrupt CO spikes superimposed. These Similar approach was applied to simulate analytical smoothing affecting other records discussed 7 in this study. It shows the baseline's enhancement due to analytical smoothing and impacting 8 other record was limited to a local maximum of 2 ppbv. 9 10 1.10. Impact of storage on CO mixing ratios 11 A 10 m long replicate section of the PLACE core was analyzed nearly two years after the main 12 analytical campaign (respectively December 2018 and February 2017) and no change in the CO 13 values was observed (Fig S14), supporting that CO production in ice during core storage is [CO] analyses on discrete ice and firn core samples (including the PLACE core) were carried out 8 at the University of Rochester (UR), with the goals of investigating the possibility of "in-extractu" 9 production of CO during melting of ice core samples from trace organics contained in the ice, as 10 well as to verify that there were no additional uncharacterized artifacts associated with the 11 continuous [CO] method. Following this, melting is started and a second expansion (Expansion 1) is performed when 1 around 0.5cm of water is present in the bottom of the vessel. This expansion occurred about 35 2 minutes after Expansion 0, with the next expansion (Expansion 2) occurring about 20 minutes 3 later, just before melting finishes, at about the time the ice bath is applied to the vessel. The last 4 three expansions (Expansions 3-5) take place about 30, 60 and 120 min after melting is complete. 5 Expansions 0-5 are reported along a time scale on Fig.3 of the manuscript. 6

Types of Tests and Data Corrections 7
To account for the procedural CO blank of the extraction line, two types of tests were performed. 8 The first is the "No-Melt" test, where the procedure is as described above except that the ice 9 sample is not melted but kept at -15°C the entire time. The average blank from all No-Melt 10 extractions was 3.2±1.7 ppbv. The second run type is the "Melt" test, involving the full procedure 11 (including melting) described above. In many cases (Table S4) The rate at which sample aliquots were expanded out of the vessel ("slow expansion" or "fast 17 expansion"; Table S4) was reduced approximately in the middle of this set of tests, as it was 18 determined that faster expansion may have allowed some water vapor to break through the 19 cryogenic trap, potentially causing an elevated response in the RGD detector (higher [CO]). The 20 results of the tests conducted using fast expansions were further corrected for this effect (1.9 21 ppbv), based on a comparison of tests run on depth-replicate firn samples using fast versus slow 22

expansions. 23
The amount of excess CO produced from melting was calculated by subtracting the starting CO 24 concentration before melting (expansion 0) from the final CO concentration after melting 25 (averaging expansions [2][3][4][5]. When the starting value (expansion 0) was not available (n=2), an 26 average of starting values (expansion 0) from similar runs was used instead. This "excess CO" 27 represents our best estimate of CO increase that is attributable solely to the melting step. 28 Finally, replicate firn samples run in "solid" and "powder" modes showed similar results (Table  29 S4), allowing to ensure that the extra sample processing (grating step) needed for ice samples 30 did not result in excess and variable CO production in the ice samples.  Table S4. Sample information for tests used to investigate the hypothesis of rapid CO production from trace organics in the ice during ice melting.

2
Firn samples were all taken from a firn core from Summit, Central Greenland (camp C14), and are denoted by a number and letter combination; the 3 number is the depth level (1)(2)(3) and the letters are the replicate sample (A-D) from that level. The "Small" Firn samples came from Firn block 3-C 4 and Firn block 3-D, which were cut in half to create "Small" Firn samples A-D. PLACE ice samples follow a similar naming scheme, with the number 5 representing the depth level and the letter indicating the replicate sample. The "PLACE ice: low CFA [CO]" samples were selected from two distinct 6 sections of the PLACE core that displayed a minimal amount of the high frequency variation in the CFA CO record, with no large CO peaks observed, 7 whereas the "PLACE ice: elevated CFA [CO]" samples were taken from a region with a large amount of the high frequency variation in the CFA CO 8 record and large CO peaks observed. Sample ice from Taylor Glacier (Antarctica) was used and is labeled by age of the ice as only a single sample 9 of each was run. Lastly, "Gas Free Ice" was produced in the laboratory and used as a control sample. The "Excess CO" reported in the table is the 1 Figure S16. Continuous CO dataset measured on two different cores collected at the NEEM site: the 2 NEEM-2011-S1 core, and the NEEM deep core (SC archive section). NEEM-2011-S1 and NEEM-SC 3 dataset overlap over the 399-409 m interval. The NEEM-2011-S1 core calibration has been revised in the 4 frame of this study. The NEEM-2011-S1 was analyzed with a Micromodule degasser. The NEEM-SC core 5 was analyzed with either a Transfer-Line or a Micromodule degasser.  The PLACE CO baseline extracted from the IGE continuous record is similar to the one based on 14 the DRI dataset: both baselines agree within their uncertainty envelopes. Replicate datasets (Fig.  15 S3) further rule out analytical drifts of the DRI or IGE CO baselines during the measurements: in 16 both laboratories, we first melted the entire core, and later some replicate sections were 17 Supprimé: 6 18 Supprimé: S16

Relationship between TOC and snow accumulation rate at Tunu13
7 Figure S26 reports the relationships between TOC and surface snow accumulation rate Tunu13.