Variability of sulfate signal in ice core records based on five replicate cores

. Current volcanic reconstructions based on ice core analysis have signiﬁcantly improved over the past few decades by incorporating multiple-core analyses with a high temporal resolution from different parts of the polar re-gions into a composite common volcanic eruption record. Regional patterns of volcanic deposition are based on composite records, built from cores taken at both poles. However, in many cases only a single record at a given site is used for these reconstructions. This assumes that transport and regional meteorological patterns are the only source of the dispersion of the volcanic products. Here we evaluate the local-scale variability of a sulfate proﬁle in a low-accumulation site (Dome C, Antarctica), in order to assess the representativeness of one core for such a reconstruction. We evaluate the variability with depth, statistical occurrence, and sulfate ﬂux deposition variability of volcanic eruptions detected in ﬁve ice cores, drilled 1 m apart from each other. Local-scale variability, essentially attributed to snow drift and surface roughness at Dome C, can lead to a non-exhaustive record of volcanic events when a single core is used as the site reference, with a bulk probability of 30 % of missing volcanic events and close to 65 % uncertainty on one volcanic ﬂux


Introduction
When a large and powerful volcanic eruption occurs, the energy of the blast is sufficient to inject megatons of material directly into the upper atmosphere (Robock, 2000).While ashes and pyroclastic materials fall rapidly to the ground because of gravity, gases remain in the atmosphere over longer timescales.Among gases, SO 2 is of particular interest due to its conversion to tiny sulfuric acid aerosols, which can potentially impact the radiative budget of the atmosphere (Rampino and Self, 1982;Timmreck, 2012).In the troposphere a combination of turbulence, cloud formation, rainout and downward transport are efficient processes that clean the atmosphere of sulfuric acid, and volcanic sulfuric acid layers rarely survive for more than a few weeks, limiting their impact on climate.The situation is different when volcanic SO 2 is injected into the stratosphere.There, the dry, cold and stratified atmosphere allows sulfuric acid layers to remain for years, slowly spreading an aerosol blanket around the globe.The tiny aerosols then act as efficient reflectors and absorbers of incoming solar radiations, significantly modifying the energy balance of the atmosphere (Kiehl and Briegleb, 1993) and the ocean (Gleckler et al., 2006;Miller et al., 2012;Ortega et al., 2015).With a lifetime of 2 to 4 years, these aerosols of sulfuric acid ultimately fall into the troposphere where they are removed within weeks.
In polar regions, the deposition of the sulfuric acid particles on pristine snow will generate an acidic snow layer, enriched in sulfate.The continuous falling of snow, the absence of melting and the ice thickness make the polar snowpack the best records of the Earth's volcanic eruptions.Hammer (1977) was the first to recognize the polar ice's propen-sity to record such volcanic history.Built on the seminal work of Hammer et al. (1977), paleovolcanism developed around this discovery and has two aims.
The first relies on the idea that the ice record can reveal past volcanic activity and, to a large extent, its impact on Earth's climate history (Robock, 2000;Timmreck, 2012).Indeed, on a millennium timescale, volcanoes and solar activity are the two main recognized natural climate forcings (Stocker et al., 2013).Based on ice records, many attempts are made to extract the climate forcing induced by a volcanic eruption (Crowley and Unterman, 2013;Gao et al., 2007Gao et al., , 2008;;Sigl et al., 2013Sigl et al., , 2014;;Zielinski, 1995).However, such an approach is inevitably prone to large uncertainty pertaining to the quality of the ice record and nonlinear effects between deposition fluxes and source emissions (Pfeffer et al., 2006).
The second aim of paleovolcanism is to provide an absolute dating scale when clear volcanic events in differently located ice cores can be unambiguously attributed to the same dated event (Severi et al., 2007).The time synchronization of different proxy records is possible, allowing the study of the phasing response of different environmental parameters to climate perturbation (Ortega et al., 2015;Sigl et al., 2015) or estimating the snow deposition over time (Parrenin et al., 2007).Whatever the aim, paleovolcanism should rely on robust and statistically relevant ice core records.
Work undertaken to date to establish a volcanic index has assumed that volcanic events are clearly identified, without any false signal from background variations induced by other sulfur sources (e.g., marine, anthropogenic).Seasonal layer counting is used whenever possible, bipolar comparison of ice sulfate records has become the method of choice to establish an absolute dated volcanic index (Langway et al., 1988).Both known and unknown events can be used to synchronize different cores.However, only a limited number of peaks, with characteristic shape or intensity, known to be associated with a dated eruption, can be used to set a reliable timescale (Parrenin et al., 2007).This restriction is partly fueled by the poor and/or unknown representativeness of most volcanic events found in ice cores.Most of the time, a single core is drilled at a given site and used for cross comparison with other sites.This approach is clearly insufficient for ambiguous events.
On a large scale, sulfate deposition is highly variable in space and mainly associated with atmospheric transport and precipitation patterns.On a local scale (ca. 1 m), variability can emerge from post-deposition processes.While sulfate is a nonvolatile species supposed to be well preserved in snow, spatial variability is induced by drifted snow and wind erosion leading to surface roughness heterogeneities (Libois et al., 2014).These effects are amplified at lowaccumulation sites where most of the deep drilling is performed (EPICA community members, 2004;Jouzel, 2013;Lorius et al., 1985).To the best of our knowledge, only one study has used multiple drillings at a given site to analyze the representativeness of the ice core record (Wolff et al., 2005).This study took advantage of the two EPICA (European Project for Ice Coring in Antarctica) cores drilled at Dome C, 10 m apart (Antarctica; 75 • 06 S, 123 • 21 E; elevation: 3220 m; mean annual temperature: −54.5 • C) (EPICA community members, 2004) to compare the dielectric profile (DEP) along the 788 m common length of the two cores.For the two replicate cores, statistical analysis showed that up to 50 % variability in the pattern of any given peak was encountered as a consequence of the spatial variability of the snow deposition.The authors concluded that ice core volcanic indices from single cores at such low-accumulation sites could not be reliable and that what was required was a network of closely spaced records.However, as mentioned in the conclusion of Wolff et al. (2005), this statistical study relied only on two records.Additionally, DEP signals are known to be less sensitive than sulfate signals for volcanic identification, and more accuracy is expected by comparing sulfate profiles.The authors thus encouraged conducting a similar study on multiple ice cores to see if the uncertainty could be reduced.
In the present study we took advantage of the drilling of five ice cores at Dome C, initially intended for the analysis of sulfur isotopes of the volcanic sulfate.Putting aside the number of records, our approach is similar in many points to the work of Wolff et al. (2005).However, it has the advantage of relying on highly resolved sulfate profiles.In addition, the spatial scale is slightly smaller as the five cores were drilled 1 m apart.The comparison of five identically processed cores is a chance to approach the representativeness of a single-core reconstruction at a low-accumulation site, the most prone to spatial variability.The representativeness of a volcanic record can be assessed by isolating the volcanic peaks in different records, as done in Wolff's work and in this study, or by a global comparison of the sulfate concentration records as proposed in Gfeller et al. (2014).In the latter case, the full individual profiles (background + the volcanic peaks) are compared to a theoretical ideal case made up of an infinite number of profiles.A similarity coefficient is then calculated between a population made up of n profiles and the infinite population.However, this approach cannot be extrapolated to discrete profiles, as in our approach, because there is a priori no ideal profile for the volcanic record.Nevertheless, the representativeness of sulfate in the Dome C record, as defined in the work of Gfeller et al. (2014), has been also calculated for comparison with this method, and the result is available in the Supplement (Fig. S1).
New constraints on the variability of sulfate deposition recorded by spatial heterogeneity at such sites are expected from the present work.Even if recent publications (Sigl et al., 2014), underline the need to use multiple records at lowaccumulation sites to overcome the spatial variability issue, such records are not always available.This lack of records adds uncertainty to the volcanic flux reconstruction based on polar depositional patterns.Our study should help to better constrain the error associated with local-scale variability and, ultimately, the statistical significance of volcanic reconstructions.The present study discusses the depth shift, occurrence of events and deposition flux variability observed in the five cores drilled.

Core drilling
The project VOLSOL (VOLcanic and SOLar natural climatic forcings), initiated in 2009, aimed at constraining the estimation of the natural part of radiative forcing, composed of both volcanic and solar contributions using ice core records of sulfate and 10 Be.In order to build a robust volcanic index including a discrimination of stratospheric events based on sulfur isotopic ratios (Baroni et al., 2008;Savarino et al., 2003), 5 × 100 m firn cores (diameter: 10 cm) were drilled in 2010/2011 along a 5 m straight line and spaced approximately 1 m apart.The drilling took place at the French-Italian station Concordia, more precisely between Concordia station and the EDC (EPICA Dome C core) drilling tent (300 m west of the EDC drilling tent).At this site, the mean annual snow accumulation rate is about 25 kg m −2 yr −1 , leading to an estimated time period covered by the cores of 2500 years.Cores were logged and bagged in the field, and temporarily stored in the underground core buffer (−50 • C) before analysis.The unusual number of ice cores drilled at the same place was driven by the amount of sulfate necessary to conduct the isotopic analysis.However, this number of replicate cores drilled 1 m apart also offers the opportunity to question the representativeness of a volcanic signal extracted from a single core per site.

Analyses and sampling
Analyses were performed directly in the field during two consecutive summer campaigns.Thirty meters were analyzed in 2011, and the rest was processed the following year.The protocol was identical for each core and the steps followed were as follows: -decontamination of the external layer by scalpel scrapping; -longitudinal cutting with a band saw of a 2 cm stick out of the most external layer; -sampling of the ice stick at a 2 cm resolution (ca.23 600 samples); -thawing the samples in 50 mL centrifuge tubes and transferring them into 15 mL centrifuge tubes positioned in an autosampler; -automatic analysis with a Metrohm IC 850 in suppressed mode (NaOH at 7 mM, suppressor H 2 SO 4 at 50 mM, Dionex AG11 column) in a fast IC configuration (2 min run) with regular calibration (every 60 samples) using a certified sulfate reference solution (Fisher brand, 1000 ppm certified).
Due to the fragility of snow cores, the first 4 m were only analyzed on a single core (Fig. 1).We will thus not discuss the variability of the Pinatubo and Agung eruptions present in these first 4 m.Concentration data are deposited in the public domain and are made freely available at the NOAA National Centers for Environmental Information, Paleoclimatology Data (https://www.ncdc.noaa.gov/data-access/paleoclimatology-data).

Peak discrimination method
As with most algorithms used for peak detection, the principle is to detect anomalous sulfate concentration peaks from a background noise (stationary or not), which could potentially indicate a volcanic event.The estimation of the background value should therefore be as accurate as possible.Using core 2 as our reference core, we observed a background average value stationary and close to 85 ppb ± 30 ppb (1σ ) at Dome C during the 2500 years of the record.However, the variability is sufficient enough to induce potential confusion on the detection of small peaks.Therefore, a stringent algorithm using the PYTHON language (accessible on demand) was developed to isolate each possible peak.The algorithm treats the full ice record by 1 m sections (ca.45-50 samples).For each meter, a mean concentration (m) and standard deviation (σ ) is calculated regardless of the presence or not of peaks in the section.Then, every value above m + 2σ is removed from the 1 m data set.A new mean and standard deviation are calculated and the same filtration is applied.Iteration runs until no more data above m + 2σ is found.At that point, m represents the background mean concentration (The resulting background estimation along core 1 is illustrated in the Supplement, Fig. S2).The process runs for each 1 m section, starting from the surface sample until the end of the core.
Then, each 1 m data set is shifted by one sample; the process is reset and the peak detection run again on each new 1 m data set.Sample shift is applied until the last sample of the first 1 m section is reached so that no bias is introduced by the sampling scheme.Every concentration data point is thus compared approximately with its 100 neighboring data points (50 on each side).Each data point isolated by the algorithm is further tested.In order to be considered as a point belonging to a potential volcanic peak, the data should be detected in a given core (i.e., the m + 2σ final threshold) in at least 50 % of the 50 runs.Additionally, the point has to be part of at least three consecutive points passing the same 50 % threshold detection.This algorithm was applied individually to each core, giving five different lists of peak.In total, 54, 51, 47, 50 and 47 peaks were detected on cores 1, 2, 3, 4 and 5, respectively.A manual detection is then required  if one wants to build a more accomplished volcanic record from several profiles, which must be based on shape criteria and not only statistical criteria.However, in the scope of this paper, no manual sorting was applied, so that the statistical assessment could rely on more objective criteria (the number of occurrences).

Core synchronization and dating
Core 1 was entirely dated with respect to the recently published volcanic ice core database (Sigl et al., 2015) using the Analyseries 2.0.8 software (http://www.lsce.ipsl.fr/Phocea/Page/index.php?id=3), and it covers the time period of 588 BCE to 2010 CE. Figure 2 shows the age-depth profile obtained for this core.A total of 13 major volcanic well-dated eruptions were used as time markers to set a timescale (bold date in Table 1).Core 1 was entirely dated through linear interpolation between these tie points.The dated core 1 was then used as a reference to synchronize the remaining four cores, using the same tie points and 10 additional peaks (nonbold date in Table 1), presenting characteristic patterns common to each core.In total, 23 points were therefore used to synchronize the cores.

Composite building from the five ice cores
Through the routine described above, the five cores are depth-synchronized using the 23 tie points, and other potential volcanic events in each core cores are detected independently.Therefore, the number of peaks detected in each core is different (between 47 and 54) and their depth (with the exception of the tie points used) is slightly different from each of the other cores due to the sampling scheme and position of the maximum concentration.After correcting the depth shift between cores, a composite profile was built by summing all Table 1.Tie points used to set the timescale and synchronize the cores.Volcanic events are named "ev x" if they are not assigned to a well-known eruption.Dating of the events is based on Sigl et al. (2015) and is given as CE dates, with negative figures referring to dates BCE. the peaks identified in the five cores.In this composite, sulfate peaks from different cores are associated with a same event as soon as their respective depth (corresponding to the maximum concentration) is included in a 20 cm depth window.This level of tolerance is consistent with the dispersion in width and shape of peaks observed (Fig. 3).A number of occurrences is then attributed to each sulfate peak, reflecting the number of times it has been detected in the five-core data set (Fig. 4).

Depth offset between cores
Depth offsets between cores are the result of the surface roughness at the time of drilling, variability in snow accumulation, heterogeneous compaction during the burying of snow layers and logging uncertainty.This aspect has been discussed previously, over a similar timescale (Wolff et al. 2005) and over a longer timescale (Barnes et al. 2006) in Dome C. Surface roughness, attributed to wind speed, temperatures and accumulation rate, is highly variable in time and space.These small features hardly contribute to the depth offset on a larger spatial scale, in which case glacial flow can control the offset between synchronized peaks, as seems to be the case at the South Pole site (Bay et al. 2010)  should be found at different depths.The offset trend fluctuates with depth, due to a variable wind speed (Barnes et al., 2006).To estimate the variability in the depth shift for identical volcanic events, we used the tie points listed in Table 1.
For each peak maximum, we evaluate the depth offset of core 1, 3, 4 and 5, with respect to core 2. To avoid logging uncertainty due to poor snow compaction in the first meters of the cores and surface roughness at the time of the drilling, we used the UE 1809 depth in core 2 (13.30m) as a depth reference horizon from which all other depth cores were anchored using the same 1809 event.For this reason, only eruptions prior to 1809 were used to evaluate the offset variability, that is 18 eruptions instead of the 23 used for the core synchronization.Figure 5  offset increases (note that the positive or negative trends are purely arbitrary and depend only on the reference used, here core 2).The maximum offset, obtained between cores 3 and 5 is about 40 cm.Such accrued offsets with depth were also observed by Wolff et al. (2005) and were attributed to the process of logging despite the stringent guidelines used during EPICA drilling.Similarly, discontinuities in the depth offset, observed by Barnes et al. (2006), were interpreted as resulting from logging errors.As no physical processes can explain a trend in the offsets, we should also admit that the accrued offset is certainly the result of the logging process.
In the field, different operators were involved, but a common procedure was used for the logging.Two successive cores extracted from the drill were reassembled on a bench to match the nonuniform drill cut and then hand-sawed meter by meter to get the most precise depth core, as neither the drill depth recorder nor the length of the drilled core section can be used for establishing the depth scale.This methodology involving different operators should have randomized systematic errors, but obviously this was not the case.Despite the systematic depth offset observed, synchronization did not pose fundamental issues as the maximum offset in rescaled profiles never exceeds the peak width (ca.20 cm) thanks to the 10 possible comparisons when a pair of cores is compared.Confusion of events or missing events are thus very limited in our analysis (see next section).Figure 5. Depth offset of 18 common and securely identified volcanic events in cores 1, 3, 4 and 5 relative to core 2. To overcome offset due to the drilling process and poor core quality in the first meters, UE 1809 (depth: ca. 13 m) is taken as the origin and horizon reference.

Variability in event occurrence
The variability in event occurrence in the five ice cores has been evaluated through the construction of a composite record (Fig. 4) and the counting of events in each core as described in the Methods section.By combining the five ice cores, we listed a total of 91 sulfate peaks (Pinatubo and Agung not included), which are not necessarily from volcanic sources.Some peaks can be due to post deposition effects affecting the background deposition or even contamination.When it comes to defining a robust volcanic index, peak detection issues emerge.The risk of misinterpreting a sulfate peak and assigning it, by mistake, to a volcanic eruption, as well as the risk of missing a volcanic peak, can be examined through a statistical analysis conducted on our five cores.We try to evaluate to what extent multiple-core comparison facilitates the identification of volcanic peaks among all sulfate peaks that can be detected in a core.To do so, we assumed that a peak is of volcanic origin as soon as it is detected in at least two cores.In other words, the probability of having two nonvolcanic peaks synchronized in two different cores is zero.It is expected that combining an increasing number of cores will increasingly reveal the real pattern of the volcanic events.All possible combinations from two-to five-core comparisons were analyzed, totalizing 26 possibilities for the entire population.The results of each comparison were averaged, giving a statistic on the average number of volcanic peaks identified per number of cores compared.The results of the statistical analysis are presented in Fig. 6.As expected, in a composite made up of one to five cores, the number of sulfate peaks identified as volcanic peaks (through being detected at least twice) increases with the number of cores combined in the composite.Thus, while only 30 peaks Figure 6.Black dots on the red line (left axis) represent the number of sulfate peaks that can be identified as volcanic peaks in a composite profile made up of n cores (with n ranging from 1 to 5).A sulfate peak appearing simultaneously in at least two cores is considered to be a volcanic peak.Blue diamonds represent the ratio of identified volcanic peaks, i.e the number of identified volcanic peaks (plotted on the left axis), relative to the total number of sulfate peaks (no discrimination criteria) in a composite made up of five cores.
In our case, the five-ice-core composite comprises 91 sulfate peaks (Agung and Pinatubo excluded).With two cores, only 33 % of them would be identified as being volcanic peaks (detected in both cores), while 68 % of them can be identified as volcanic events using five cores.
can be identified as volcanic from a two-core study, a study based on five cores can yields 62 such peaks.The five-core comparison results in the composite profile given in Fig. 4a.
The initial composite of 93 peaks is reduced to 64 volcanic peaks (Pinatubo and Agung included) after removing the single peaks (Fig. 4b).Each characteristic of the retained peaks is given in Table 2.The main conclusion observing the final composite record is that only 17 of the 64 peaks were detected in all of the five cores and 68 % of all peaks were at least present in two cores.On the other end of the spectrum, two-core analysis reveals that only 33 % (30 peaks on average) of the peaks are identified as possible eruptions.A twocore comparison still presents a high risk of not extracting the most robust volcanic profile at low-accumulation sites, a conclusion similar to that of Wolff et al. (2005).Surprisingly, it may also be noted that this five-core comparison does not result in an asymptotic ratio of identified volcanic peaks, suggesting that five cores are not sufficient either to produce a full picture.High-accumulation sites should be prone to less uncertainty; however, this conclusion remains an a priori that still requires confirmation.Large and small events are not equally affected by these statistics.Figure 7 shows that the probability of presence is highly dependent on peak flux and the risk of missing a small peak (maximum flux in the window (f + 2σ : f + 5σ ), f be- ing the background average flux) is much higher than the risk of missing a large one (maximum flux above f + 8σ ).However, it is worth noting that major eruptions can also be missing from the record, as has already been observed in other studies (Castellano et al., 2005;Delmas et al., 1992).The most obvious example in our case is the Tambora peak (1815 AD), absent in two of our five drillings while presenting an intermediate to strong signal in the others (Fig. 8).The reason for the variability in event occurrence has been already discussed by Castellano et al. (2005).In the present case of close drillings, long-range transport and large-scale meteorological conditions can be disregarded due to the small spatial scale Table 2. Sulfate peaks (maximum concentration in nanograms per gram and flux of volcanic sulfate deposited in kilograms per square kilometer) considered to be volcanic eruptions based on the statistical analysis of the five cores.Flux is calculated by integrating the peak, using the density profile obtained during the logging process.Volcanic flux values are corrected from background sulfate (calculated separately for each sulfate peak).Zero stands for non-detected events in the cores.Agung (3.77 m) and Pinatubo (1.52 m) were not included in the statistical analysis because they were analyzed only in core 1 and thus are marked as not applicable (n/a).The estimation of the average volcanic flux takes into account undetected peaks, for which the flux is considered to be 0. The relative error in the flux (estimated as 10 %) takes into account the IC measurement relative standard deviation (below 4 % based on standards runs), the error in firn density (relative error estimated as 2 %) and the error in sample time length (10 %).The last column displays data obtained from Castellano et al. (2005) for identical volcanic peaks.For similar peaks Castellano's flux generally falls into the average flux +40 % uncertainty, sometimes exceeding this value.Dates are given as CE dates, with negative numbers indicating dates BCE.   of our study; the snow drift and surface roughness is certainly the main reasons for missing peaks.The fact that two events as close to one another as UE 1809 and Tambora are so differently recorded indicates that post-depositional effects can affect the recording of eruptions very variably in time and space.

Variability in signal strength
To compare peak height variability, detected peaks were corrected by subtracting the background from peak maxima.We considered C i /C mean variations, C i being the SO 2− 4 maximum concentration in core i (1 to 5) and C mean being the mean of these concentrations for the event i. C i is considered zero if the peak is not detected in a core.For concentration values, positive by definition, the log-normal distribution is more appropriate; geometric means and geometric standard deviations were used, as described by Wolff et al. (2005) (Table 3).In our calculation, the geometric standard deviation based on five cores is 1.49; in other words, the maximum concentration of a peak in one core is uncertain by 49 %.This factor is completely in agreement with the one obtained in Wolff et al. (2005) (1.5).Having n cores allows for a reduction in the uncertainty on the mean (standard error of the mean) by a factor 1/ √ n.The peak height mean obtained from five cores is therefore uncertain by 22 %.Comparing peak maxima induces a bias related to the sampling method: with a 2 cm resolution on average, a peak's height is directly impacted by the cutting, which tends to smooth the maxima.Comparing the total sulfate deposited during the event is more appropriate.Pursuing a similar approach but reasoning on the basis of mass of deposited sulfate rather than maximum concentration (and considering F i /F mean , F i being the mass flux of peak i), the obtained variability is higher than previously.The uncertainty in the flux for one measurement is 65 % (based on the standard deviation of the mean), and the uncertainty of the mean (standard error of the mean) is therefore close to 30 %.The difference in the signal dispersion between the two approaches rests on the fact that the peak maximum has a tendency to smooth the concentration profile as a consequence of the sampling strategy.This artifact is suppressed when the total mass deposited is considered.

Conclusion
This study confirms in many ways previous work on multiple drilling variability (Wolff et al., 2005).As already discussed, peaks flux uncertainty can be significantly reduced (65 to 29 %) by averaging five ice core signals.A five-core composite profile was built using the criterion that a peak is considered volcanic if present in at least two cores.We observed that the number of volcanic peaks listed in a composite profile increases with the number of cores considered.With two cores, only 33 % of the peaks present in the composite profile are tagged as volcanoes.This percentage increases to 68 % with five cores.However, we did not observe an asymptotic value, even with five cores.A single record at a low-accumulation site is therefore very unlikely to be a robust volcanic record.Of course, peaks presenting the largest flux are more likely to be detected in any drilling, but the example of the Tambora eruption shows that surface topography is variable enough to erase even the most significant signal, although this occurs rarely.This variability in snow surface is evidenced in the depth offset between two cores drilled less than 5 m from each other, as peaks can easily be situated 40 cm apart.
At low-accumulation sites such as Dome C, where surface roughness can be on the order of the snow accumulation and highly variable, indices based on chemical records should be considered with respect to the timescale of the proxy studied.Large timescale trends are only a little sensitive to this effect.By contrast, a study on episodic events such as volcanic eruptions or biomass burning, with a deposition time on the order of magnitude of the surface variability scale should be based on a multiple-drilling analysis.A network of several cores is needed to obtain a representative record, at least in terms of recorded events.However, although lowered by the number of cores, the flux remains highly variable, and the mean flux obtained from five cores is still uncertain by almost 30 %.This point is particularly critical in volcanic reconstructions that rely on the deposited flux to estimate the mass of aerosols loaded in the stratosphere and, to a larger extent, the climatic forcing induced.Recent reconstructions largely take into account flux variability associated with a regional pattern of deposition, but this study underlines the necessity of not neglecting local-scale variability at low-accumulation sites.Less variability is expected with a higher accumulation rate, but this still has to be demonstrated.Sulfate flux is clearly one of the indicators of eruption strength, but due to transport, deposition and post-deposition effects, such a direct link should not be taken for granted.
With such statistical analysis performed systematically at other sites, we should be able to reveal even the smallest volcanoes imprinted in ice cores, extending the absolute ice core dating, the teleconnection between climate and volcanic events and improving the time resolution of the mass balance calculation of ice sheets.

Figure 1 .Figure 2 .
Figure 1.Sulfate profiles on the five replicate cores obtained during a drilling operation at Dome C -Antarctica in 2011.Data are available at https://www.ncdc.noaa.gov/data-access/paleoclimatology-data.

Figure 3 .
Figure 3. Kuwae (a), Krakatua (b) and Tambora (c) sulfate concentration profiles after depth synchronization.All peaks are within a 20 cm uncertainty, enabling us to clearly attribute each occurrence to a single event.

Figure 4 .
Figure 4. Panel (a): composite sulfate peak profile deduced from our statistical analysis of the five cores using our detection peak and synchronization algorithms (see text).The numbers indicate the number of time a common peak is found in the cores.Unnumbered peaks: peaks found only in single core.Panel (b): same as (a) without the single detected peaks.All the remaining peaks are considered to be volcanic eruptions.See Table2for details.

Figure 7 .Figure 8 .
Figure 7. Peaks probability to be detected in two, three, four or five cores, as a function of their flux.The three categories of fluxes are defined by peak flux value, relatively to the average background flux, and quantified by x time (2, 5 and 8) the flux standard deviation (calculated for a 30 ppb standard deviation in concentrations).At flux above background flux + 8σ , there is a 90 % chance of detecting the volcanic peak in each core of a population of five cores.On the other hand, at flux below background flux + 5σ , there is a 60 % probability of the volcanic peak being detected in two cores only among the five-core population.This highlights that replicate cores are particularly useful to avoid missing small to intermediate peaks in a record.

Table 3 .
Statistics on sulfate signal for identical peaks in cores 1, 2, 3, 4 and 5. Geometric standard deviations are calculated on peak heights (i.e., maximum concentration reached in nanograms per gram) and on peak sulfate flux (i.e., total mass of volcanic sulfate deposited after the eruption).Background corrections are based on background values calculated separately for each volcanic event.