Interactive comment on “ Exploitation of chemical profiles by conjugate variable analysis : application to the dating of a tropical ice core ( Nevado Illimani , Bolivia ) ”

Introduction Conclusions References


Conclusions References
Tables Figures

Back Close
Full

Abstract
Ice core dating is a key parameter for the interpretation of the ice archives.However, the relationship between ice depth and age can generally not be easily established and requires to combine a large number of investigations and/or modeling effort.This paper presents a new approach of ice core dating based on conjugate variable (depth and spatial frequency) analysis of chemical profiles.The relationship between the depth of a given ice layer and the date it was deposited is determined using ion concentration depth profiles obtained along a one hundred-meters deep ice core recovered in the summit area of the Nevado Illimani (6350 m a.s.l.), located in the Eastern Bolivian Andes ( 16• 37 S, 67 • 46 W).The results of Fourier conjugate analysis and wavelet tranforms are first compared.Both methods are applied to nitrate concentration depth profile.The resulting chronologies are checked by comparison with the multi-proxy year-by-year dating published by de Angelis et al. (2003) and with volcanic tie points, demonstrating the efficiency of Fourier conjugate analysis when tracking the natural variability of chemical proxies.The Fourier conjugate analysis is then applied to concentration depth profiles of seven other ions thus providing information on the suitability of each of them for dating studies of tropical Andean ice cores.

Introduction
Ice core dating, i.e. establishing an univocal relationship between ice depth and age, is not easy, although fundamental for the interpretation of ice archives.Current methods generally combine various approaches like visual stratigraphy, seasonal features allowing to identify annual snow layers, location of specific physico-chemical patterns previously identified in well dated sedimentary or ice archives, and ice flow or inverse modeling.
Annual layer counting is well suited to the study of high frequency atmospheric processes but requires a certain amount of refinement, due to uncertainty in identifying in-Figures

Back Close
Full dividual layers.In particular, annual signal may be distorded or smoothed by changes in snow accumulation rate and processes involved in snow deposition and firn aging: snow blowing, wind pumping, water vapor or gradient diffusion and interaction between reactive species may significantly modify the initial registers.Dating uncertainty increases with depth, making it very useful to identify, if any, paleo-events allowing to synchronize different ice cores and potentially to check the current dating for time periods longer than one or two centuries (Narcisi et al., 2005;Traufetter et al., 2004).The very large number of studies undertaken over polar ice caps since the end of the 1950's has helped to considerably improve our understanding of processes leading to the building of polar ice records and their dating at least for central areas of Antarctica and Greenland.The comparison between deep ice cores recovered in Central Greenland and in Antarctica as well as marine sediment ice cores has provided even more accurate ice chronologies for time periods covering the most recent decades or the past few centuries as well as several hundred of kyr (Ruth et al., 2007;Lemieux et al., 2010).A typical example is the multi-parameter approach combining visual stratigraphy, identification of seasonal maximum along chemical profiles, and age matching with specific climatic and volcanic events used as tie points which provided a year-by-year dating of the GISP2 ice core for the past 100 kyr (Meese et al., 1997).In addition to modeled time scales, new Greenland time scale has been more recently developped back to 60 kyr before 2000 AD, based on visible stratigraphy and the identification of annual layers along a large set of high resolution records of ice chemical and physical properties (Andersen et al., 2006;Svensson et al., 2005;Rasmussen et al., 2008;Svensson et al., 2008).The thorough comparison of annual layer counting carried out by independented experienced investigators greatly improved the dating accuracy but is a tedious time-consuming effort that highlights the need for automated and less subjective method.
The study of Andean ice records started in the mid 1980s with paleo-climate reconstructions based on the comparison of water stable isotopes profiles with Antarctic or Greenland data (Thompson et al., 1998(Thompson et al., , 2003)).More recently, studies have focused Figures

Back Close
Full specifically on understanding the atmospheric natural variability at the annual and interannual time scale (de Angelis et al., 2003;Correia et al., 2003;Vimeux et al., 2008).However, logistical constraints and harsh weather conditions prevailing over the high altitude cold glaciers of the Andean Cordillera make it often difficult to get the full set of glaciological and meteorological data required by the understanding of processes leading to the building of Andean ice archives.Thus, improving the dating of the cores by developing reliable dating tools remains a key challenge.Particular attention must be paid to sublimation processes and snow remobilization by wind which is suspected to play a prominent role at sites influenced by major cyclonic activity or seasonal convection.Moreover, the rich and intricate fingerprint of Andean volcanism superimposes the long range influence of cataclysmic eruptions commonly used as time markers along polar ice cores.Nevertheless, most Andean ice cores were recovered in the tropical belt where the precipitation system is characterized by the yearly alternation of a warm wet season (around NDJ) and a colder dry season (around JJA) accompanied by significant changes in agricultural practices, soil emissions, and air mass trajectories.Biomass burning, increased aridity, and soil dust mobilization send large amounts of aerosol and gases into the tropical atmosphere during the dry season at the same time as limited water vapor condensation and/or riming occur at the surface of falling snow crystals.As a consequence, the concentration of chemical species in fresh snow should be significantly higher during the dry season than during the wet season.The increasing effect of dry deposition or snow sublimation on the concentration of irreversible species in aging surface snow are also likely to be more pronounced during the dry season which justifies taking into account seasonal patterns as potential dating tools along ice cores extracted from tropical Andean cold glaciers.However, seasonal trends in ice cores may be overlapped by other components related to the inter-annual variability of source properties (location, emission intensity), of large scale meteorological patterns, and of precipitation occurrence at these high altitude sites, so that an interesting step was to make the year-by-year counting independent of the feeling of the observer and to minimize the risk of cumulative errors.Introduction

Conclusions References
Tables Figures

Back Close
Full Various attempts based on spectral analysis have been developed for automatically exploring sedimentary archives.Preto et al. (2004) for instance used Time-Frequency analysis to decipher the influence of periodic features in Earth current orbit rotational configuration (precession index) on the building of sedimentary archives.made of well defined lithofacies unambiguously identified by field observation and additional laboratory studies.Smith et al. (2009) present an automated trace element laminae peak counting procedure based on spectral and wavelet analysis and providing reliable and accurate dating of speleothems with rather regular grow rates.Wheatley et al. (2012) developed a dating method using the autocorrelation function for large core subsection taken at different depths and where missing values or unclear annual signals are reconstructed assuming continuous annual layer thickness.An annual layer detection routine based on the framework of hidden Markov models has been applied by Winstrup et al. (2012) to a stable water isotope record and visual stratigraphy of two Greenland cores showing well preserved annual signals at great depth and where manual counting has already been carried out.Once all parameters included in the layer model have been optimally chosen by the algorithm according with the layer signal in the data (assumed to be well constrained), the algorithm reproduces pretty well the manually counted timescales.
All these studies focuse on archives containing well marked time related sequences.In particular, automated approaches for annual layer counting along ice cores have been applied to records showing strong and well constrained seasonal variability and gained at sites fully investigated like Antarctic Peninsula or Central Greenland.The context of Andean ice cores is quite different: dating tools remain very limited and their accuracy is hard to check.
This paper presents a new approach to ice core dating based on the time-frequency analysis of chemical depth profiles.The method was checked along the upper 100 m of a 137 m deep ice core recovered in 1999 close to the summit of the Nevado Illimani in the Eastern Bolivian Andes (6350 m a.s.l., 16 • 37 S, 67 • 46 W).The Illimani ice core was chosen to check the validity of our approach because: (1) continuous high reso-Introduction

Conclusions References
Tables Figures

Back Close
Full   et al., 2003).This dating made it necessary to combine the seasonal signals of three wind-borne dust markers (mass of insoluble particles, calcium, and aluminium, Correia et al., 2003) and two nitrogen ionic species (nitrate and ammonium).Two pure sulfuric acid peaks were observed at 9.13 m and 33.43 m.They were attributed to the stratospheric fallout in early 1993 and 1964, respectively, of the Pinatubo (June 1991) and Mont Agung (March 1963) eruptions.Debris corresponding to the arrival of tropospheric plumes during the 1991 and 1963 dry seasons were found at 10.8 and 33.4-33.8m, respectively, and used as check points.Dating uncertainty was estimated to be in the range of ±1 yr down to 41 m (1947 AD) and ±2 yr at the bottom of the core section (1919 AD).Data were processed in two ways: in a first step, the time-depth relationship was deduced from the time frequency analysis of the nitrate profile based on Fourier conjugate analysis and implemented to construct the ice chronology.Then a similar approach was done using wavelet transforms, traditionally used for climatic paleo-records once core chronology is achieved.We first present the raw data and the method we used to compensate for signal bias induced by the experimental sampling rate which was too low below 68 m and we propose and discuss a preliminary dating.The chronologies resulting from the two spectral methods are compared with the chronology proposed for the upper part of the core by de Angelis et al. (2003) and the volcanic tie points defined by these authors.This provides an assessment of the reliability of the two spectral approaches.In the second step, the Fourier conjugate transform is applied to the concentration depth profiles of seven other ions providing seven other independent Introduction

Conclusions References
Tables Figures

Back Close
Full ice chronologies.The comparison of these 8 independent dating results between each other and with the tie points used by de Angelis et al. (2003) is discussed in terms of method and marker suitability.

Raw data
In a first step, the method is described and applied to the nitrate concentration depth profile along the Illimani core, c NO − 3 (n).Nitrate is one of the major ions detected in the core and is less sensitive to post-deposition effects than other biomass markers like formate.Moreover, it is expected to have a strong annual component and to be not markedly influenced by sporadic events like volcanic eruptions or dust storm which can disrupt the seasonal pattern of other ions.Therefore it was considered as most suitable for checking our method.Raw data of c NO − 3 (n) are presented in Fig. 1.To counterbalance firn densification, all depths and layer thickness are expressed in meters of water equivalent (m w.e.).The depth d (m w.e.) expressed in meters of water equivalents and corresponding to a given depth d (m), measured along the core and expressed in meters, is given by: where δ(i ) is the density of the firn depth increment ∆d (i ).
In the left part of Fig. 1, we have reported nitrate concentrations as a function of depth.The three panels in the left part of the figure correspond to selected zoomed sections at 25-30 m w.e., 42-47 m w.e., and 62-67 m w.e. from top to bottom, respectively.The nitrate depth profile clearly shows an alternation of relative maxima (dry seasons) and minima (wet seasons) corresponding to its annual cycle.The spatial frequency 1 of this cycle increases from the surface to 48 m w.e.This increase of the Introduction

Conclusions References
Tables Figures

Back Close
Full The apparent decreasing trend of spatial frequency below 48 m w.e. is due to a change in sampling resolution resulting in inadequate sample length.This undersampling problem and the approach to address it are discussed in Sect.4.2.5.

Principle
As stated above, our investigation is based on the hypothesis that concentrations of chemical species entrapped in falling snow are periodic in time and that the period is known.We assume that the variability of the chemical composition of tropical Andean precipitation is controlled by strong seasonal changes in atmospheric properties with a time period of one year.However, and in spite of a good preservation of aging snow, the time dependency of chemical signals cannot be directly inferred from concentration depth profiles due to the interannual variability of annual net accumulation rate and to ice thinning through firn densification and ice flow: the annual layer thickness, a(d ), varies with depth.A relationship can be write between a(d ) and the spatial frequency f s (d ): Introduction

Conclusions References
Tables Figures

Back Close
Full The evolution of the spatial frequency associated with the variable periodicity as a function of depth is determined using Fourier conjugate variable analysis 2 Flandrin  (1993) and orthogonal Morlet wavelet transforms Goupillaud et al. (1985), applied to the nitrate concentration depth profile.The spatial frequency f s (d ) establishes the relationship linking depth and time delay, t(d ), where t(d ) is the difference between the drilling date and the deposition date of the layer located at the depth d .
In this way, we get the time-depth relationship along the part of the Illimani ice core where continuous ion concentration profiles are available, i.e. from the surface to 76.4 m w.e.depth.A comparison is made between results derived from the Fourier conjugate variable analysis, the wavelet transforms and a previously published yearby-year dating approach to evaluate the accuracy of spectral analysis methods applied to chemical species.Independent chronologies are then established using nitrate, formate, oxalate (indicators related to continental biomass), calcium, magnesium, potassium, and sodium (erodible soils) and sulfate (more influenced by volcanism).

Determination of the time-depth relationship
Let t(d ) be the relation linking the time t with the depth d .This function is developed at the first order in the neighborhood of the depth d , assuming that for small variations3 , the terms at the higher order are insignificant.
Full Let a T (d ) be the layer thickness corresponding to the depth increase, during a period of T years at the depth d .a T (d ) is given by4 : To this spatial period, we relate the spatial frequency using: Let d min be the depth origin of the determination of the spatial frequency and t(d ) the time delay from the drilling date to the depth deposition time.Using Eq. ( 5) we get: giving the relation between time and depth from the spatial frequency.This relation is the central point of our method: it shows that the time depth relationship can be determined from the measurement of the spatial frequency as a function of depth.
In Andean sites, the net accumulation rate corresponds to a periodicity of T = 1 yr.So the relationship between time and depth is given by Eq. ( 6) where T = 1.

Data processing
As stated above, data processing is first applied to the nitrate concentration depth profile.Introduction

Conclusions References
Tables Figures

Back Close
Full (n) is given for a set of depth d not regularly distributed from 0 to 76.4 m w.e.This concentration depth profile has high amplitude spikes related to particularly intense biomass burning events during the dry season and its mean value is different from zero.
Conjugate variables analysis requires samples regularly distributed with depth, without high impulses and with a mean value equal to zero.In order to get a concentration profile suitable for conjugate variables analysis the data set is oversampled by a spline interpolation on 2048 samples regularly distributed in depth.As shown in Fig. 3 where c NO − 3 (n) raw data are reported along with spline interpolation, the interpolation fits pritty well the original data.Then, the highest pulses are removed by thresholding.The threshold is determined using the cumulative distribution function of measured concentrations.Optimizing the thresholding effect on measured concentrations led us to fit the threshold to the value corresponding to the 0.98 level of the cumulative distribution function.Finally, we remove the mean concentration value.
After spline re-sampling, thresholding, and mean value removal, we estimate the averaged concentration spectrum (Fig. 4).The low frequency high level visible in Fig. 4 is due to the red character of geophysical noise.These dominant low frequencies would perturb the tracking of the spatial frequency by conjugate variables analysis.They are eliminated by applying a high-pass filter of cut-off frequency f b (Fig. 4).The spatial low frequency estimate is sensitive to the value of this cut-off frequency.Section 3.2.4explains how the cut-off frequency is optimized using volcanic markers.
The normalized new concentration obtained after this pre-processing is called S X (l ), where X is the component under investigation and l is the sample index (l varying from 1 to 2048).The initial concentration depth profile, C NO − The sampling period P s resulting from the spline interpolation is: This sampling period enables us to represent the concentration inside the spatial frequency band from 0 to the "Shannon spatial frequency", f Shannon :

Estimate of the spatial frequency
The spatial frequency is obtained by conjugate variable analysis, which provides the energy distribution of the concentration signal in the depth-spatial frequency plane.The watershed of this depth-spatial frequency representation (Flandrin, 1993;Cohen, 1994) gives the spatial frequency f s (d ) as a function of depth d .There are many DFR.Representations based on Fourier transforms provide spectrograms, while wavelet transforms provide scalograms (Daubechies, 1990;Rioul and Duhamel, 1992).In the next section, we present spatial frequency estimates deduced first from the spectrogram and then from the scalogram.

Spectrogram and scalogram
The spectrogram is obtained by pre-windowing the signal around a given depth, calculating its Fourier transform, taking the square modulus and doing the same for every depth.Applied to our signals, the spectrogram gives a DFR of the signal energy spectrum in the depth-spatial frequency domain.The DFR spectrogram is given by:

Conclusions References
Tables Figures

Back Close
Full Two orthogonal Morlet wavelets (Goupillaud et al., 1985) with S successive compression factors (scales) are used to calculate the scalogram.At each depth and at each scale we project the concentration S NO − 3 on the vector of the S orthogonal Morlet wavelets.The power at each scale is the sum of the square of the projection on the two orthogonal Morlet wavelets.This gives the power of the concentration signal on the depth-scale plane (depth-scale representation).From the relation between the scales and the spatial frequencies we get, by a spline transform, the DFR given by the power of the scalogram in the depth-spatial frequencies plane called SCAL X (d , f s ).
Figure 6 shows the spectrogram and the scalogram of nitrate concentration in the depth spatial frequency plane.The spectrogram is calculated for N consecutive samples.Based on successive tests, we determine that the best plots of the spectrogram are obtained with N = 128 samples corresponding to a layer thickness of 4.775 m we.
For the scalogram we use S = 35 scales that covers continuously the spatial frequency band of the concentration signal.
The depth trend of energy main value, and thus of spatial frequency, is tracked by the green "streak" on the spectrogram and on the scalogram shown in Fig. 6 (bottom and top, respectively).Spatial frequency tends to increase from approximately Introduction

Conclusions References
Tables Figures

Back Close
Full in Fig. 7.The sampling spatial period of this simulated concentration is 0.0625 m w.e.corresponding to the Shannon spatial frequency: From 0 to 48 m w.e., the frequency of the simulated signal increases from 2 to 8 (m w.e.) −1 .Beyond 48 m w.e., the signal frequency continues to increase but due to the aliasing, it appears to decrease.This simulation generates the same trend as the one observed for S NO − 3 in Fig. 6 and corroborates our interpretation of the decreasing streak in the S NO − 3 spectrogram and scalogram as spectrum replication occuring very likely when the signal frequency becomes higher than the Shannon frequency.
Nevertheless a question remains.We have indicated that concentrations used after the spline resampling have a Shannon spatial frequency of 13.4 (m w.e.) −1 while spectrum replication takes place at 8.59 (m w.e.) −1 (see Sect. 4.2.5).What can explain such a discrepancy?It is explained by the sampling procedure used for chemical analyses.The initial core sampling was continuous: the core was cut into 1697 successive ice sections and analyses were performed on every ice section once melted.This analytical procedure is equivalent to an averaging sampler.If we assume that the spectrum replication takes place at 8.59 (m we) −1 it implies that ice sections used for individual measurements below 48 m w.e. were of 5.84 cm w.e.thick which compares pretty well with the true sample thickness (5.3 ± 1.5 cm w.e.) at this depth.As the sampling period of the spline re-sampling is lower than the sampling period imposed by experimental conditions, the aliasing appears at the Shannon spatial frequency imposed by the chemical analysis.Introduction

Conclusions References
Tables Figures

Back Close
Full where SPECT X (d , f s ) and SCAL X (d , f s ) are the spectrogram and the scalogram corresponding to the component X (here nitrate), and are expressed as functions of depth and of spatial frequency.
Figure 8 shows spatial frequencies estimated from the spectrogram (left) and the scalogram (right).Estimated spatial frequencies increase from 0 to 48 m w.e. and then decrease due to spectrum replication.Discontinuities are visible on both panels, however their number is significantly greater on the scalogram estimate than on the spectrogram.These discontinuities come from transients.This is because the concentration is made of successive events of short duration and varying amplitude and their occurence is modulated by a seasonal periodicity.It is well known that the wavelets are very sensitive to transients (Mallat et al., 1992), and, as discussed below, the greater number of discontinuities observed in the scalogram led us to consider the spectrogram as a better tool than the scalogram for tracking the spatial frequency.

Discontinuity correction with a median filter
Median filters are used to eliminate the discontinuities in spatial frequency estimates, f s (d ) (Coquerez and Philipp, 1995).We use here a median filter of length M = 2m + 1.The output of median filters applied to the spatial frequency estimates, f s (d ), provided by the spectrogram and the scalogram is given by:

Conclusions References
Tables Figures

Back Close
Full The corresponding results are shown in Fig. 9.

Filters optimization
In the pre-processing step (see Sect. 4.1), we have introduced a high-pass filter with a cut-off frequency f b .Here we use a median filter of length M. To determine the optimum values of the high pass filter cut-off frequency, f b , and the median filter length, M, we use two volcanic references: the depths where the volcanic debris attributed to tropospheric deposits from the Pinatubo and the Mount Agung eruptions were found (de Angelis et al., 2003).The Tambora eruption is not taken into account here because the corresponding layer is located below the aliasing depth.Eruption dates were calculated for the following set of f b and M values: The optimum values of f b and M are the ones that minimize the distance between estimated dates and historical eruption dates using the criterion:

De-aliasing
As stated above, the decreasing trend of f s (d ) in Figs.date corresponding to a given depth d are given by: and the corresponding dates by: The chronologies respectively derived from the nitrate spectrogram and scalogram are shown in Fig. 11.The drilling date was fixed at 1999.39.Both chronologies are globally coherent but local discrepancies are visible, particularly near the aliasing edge.This raises the key question of the chronology reliability.

Evaluation of the chronology inferred from NO − 3 concentration depth profile
Evaluating the chronology reliability (accuracy) is a key step to validate our method.For that, we propose three approaches.The first approach refers only to our chronology.In the second approach we use the volcanic markers.Finally, we compare our chronology with the independent year-by-year chronology used by de Angelis et al. (2003),and Correia et al. (2003) for the upper part of the core.

Direct evaluation of our method
With the spectrogram and the scalogram, we have determined the time delay estimate t(d ).From the nitrate concentration depth profile S NO − The spectrograms of S NO − 3 ( t) corresponding to Fourier conjugate analysis (spectrogram, top panel) and wavelet transforms (scalogram, bottom panel) are presented in Fig. 12.
These spectrograms show two very different parts: before and after 160 yr (before the drilling date, i.e. 1999 AD).
The part of the spectrograms corresponding to time values greater than 160 yr are associated to the aliased part of the data and, due to the aliasing, all energies are concentrated on frequencies lower than 1 yr −1 .Thus, we cannot use this part of the spectrograms to check the quality of our method.
For time values lower than 160 yr, i.e., the non-aliased part of the data, the spectrogram of S NO − 3 ( t) obtained from the dating with the spectrogram (top panel) clearly shows a track centered around 1 yr −1 .This indicates that the results of the Fourier conjugate analysis are in agreement with our working hypothesis that aimed to find a seasonal triggering in the nitrate concentration depth profile.On the opposite, an increasing trend from 0.8 yr −1 to 1.2 yr −1 is visible between the origin and 160 yr on the spectrogram of S NO − 3 ( t) obtained from the dating with the scalogram (bottom panel).This strongly suggests that wavelet transforms cannot be used to unambiguously detect seasonal variations.

Evaluation of the chronologies deduced from NO − 3 depth profile by comparison with volcanic tie points
The historical dates of Pinatubo, Mt Agung and Tambora tropospheric plume arrival and the date deduced from the use of the spectrogram (Date 1) and the scalogram (Date 2) are given in Table 1.As a general rule, specific attention must be paid to the use of volcanic eruptions as tie points for the dating of Andean ice cores.Major volcanic eruptions can be clearly identified in Antarctic ice cores through the arrival of large amounts of sulfuric acid from the stratosphere.Nevertheless, the formal identification of a given individual event requires specific micro-physical and micro-chemical inves-Introduction

Conclusions References
Tables Figures

Back Close
Full tigations of ashes to be confirmed unless detailed matching with another well dated ice core can be provided.In tropical South America, the quiescent or eruptive activity of Andean volcanoes competes with the influence of eruptions taking place farther north in the tropical belt.Combined with the lack of information on dispersion and removal parameters as well as on the relative importance of tropospheric transport and stratosphere-troposphere exchanges, this makes uncertain the use of general volcanic markers alone to relate a given peak to a given eruption.The depth corresponding to the plume arrival attributed to Mount Agung eruption is in agreement with 37 Cs and tritium peaks at 33.5 m (21.7 m w.e.) associated with stratospheric fallout in 1964 of major nuclear bomb tests performed in the Northern Hemisphere.
The dates obtained for the Pinatubo and Mt Agung eruptions with the two approaches are very close to the eruption historical dates.Note that these references are also used to optimize the high pass filter cut-off frequency and the median filter length.On the other hand, the dates deduced from the spectrogram and the scalogram for the arrival of volcanic debris attributed to the Tambora eruption by de Angelis et al. (2003) are significantly different from the historical date, i.e. 1815.These may be due to the fact that this event is located in the aliased part of the data, but they must also be considered keeping in mind chronology discrepancies between previous studies related to the choice of time references considered to be "certain".The chronology of a nearby ice core based on the counting of pseudo-seasonal layers and the location of major peaks in the electrical conductivity profile Knüsel et al. (2003) locates at 99.9 m w.e. the arrival of a major unknown eruption dated in 1258 AD in Antarctica.The choice of this volcanic anchor layer appeared to be congruent with visual counting.However this chronology was largely reassessed by Kellerhals et al. (2010).Considering thallium as an unambiguous volcanic marker, these authors used the location of major thallium peaks combined with ice modelling and preliminary (14) C data to locate the 1258 AD horizon at 90.7 m w.e., i.e. 9.2 m w.e.higher, which implies a gap of more than 300 yr between both chronologies and underlines the potential bias induced by subjective factors in visual stratigraphy and the choice of tie points.Introduction

Conclusions References
Tables Figures

Back Close
Full

Comparison with a year-by-year multi-proxy counting approach
We compare the time scales respectively deduced from the S NO − 3 spectrogram and scalogram for the upper part of the core (0-36.4m w.e.) with the chronology inferred from the multi-proxy year-by-year dating approach proposed by de Angelis et al. (2003). The

Conclusion on the application of Fourier and wavelet transforms to the NO − 3 concentration depth profile
In the previous sections we presented: (i) a direct evaluation of our chronology; (ii) the comparison of the dates deduced from spectrogram and scalogram analysis for volcanic layers attributed to well known cataclysmic eruptions and the true eruption dates; (iii) the comparison of time scales provided by scalogram and spectrogram analysis with a previously published chronology based on a multi-proxy year-by-year counting approach.This allows us to conclude that our method provides chronologies consistent with available information.However, the seasonal variability of nitrate concentrations is clearly identified only when using the Fourier conjugate analysis.Moreover, the chronology deduced from the spectrogram shows a better general agreement with the published by de Angelis et al. (2003) time scale than the chronology inferred from the scalogram.As shown in Fig. 8, this may be explained by the much greater rate of dis-

Conclusions References
Tables Figures

Back Close
Full continuities introduced when estimating spatial frequencies from the scalogram.It is not relevant to directly evaluate the quality of our method in the aliased part of the data.However, it must be noted that although the Tambora is inacurrately dated, we were able to detect spectrum replication due to the sampling procedure, to reverse the aliasing and to propose a preliminary time scale.Although not accurate if compared with what may be expected from high resolution concentration depth profiles, this time scale is also more reliable than the year-by-year counting performed on the under-sampled part of the core as proposed by Ramirez et al. (2003).Indeed, these authors tracked seasonal variations down to 66 m w.e. but with a mean sample length of 5.8 cm, i.e. still greater than in this study that should have led to significantly overestimate accumulation rates (see also Fig. 1).Using the sulfate layer attributed to the Tambora eruption as an absolute tie point, they estimated the uncertainty of their multi-proxy approach to be ±20 yr at 66 m w.e.This is in the range of the date difference found here but fails to take into account the additional error introduced by the too low sampling rate which is likely to reach several tens of years at this depth.
In the next section, the spectrogram is applied to concentration depth profiles of seven other ions, making it possible to compare and discuss the suitability of each ion as a dating proxy.

Suitability of chemical markers for Illimani ice core dating
The seven ions studied here and compared with nitrate (NO − 3 ), are sulfate (SO 2− 4 ), sodium (Na + ), calcium (Ca 2+ ), magnesium (Mg 2+ ), potassium (K + ), oxalate (Ox), and formate (For).They were chosen for the reservoirs they are expected to originate from and for their emission seasonality.The soluble part of continental input at the Illimani site is mainly regional and consists of gypsum (CaSO 4 ) and halite (NaCl).Strong acidic gases of volcanic origin react with wind-blown minerals and volcanic ashes leading to the formation of secondary gypsum, halide salts and lower amounts of other sodium, magnesium and potassium salts.Halite may also originate from the large CPD 9, 3399-3447, 2013

Exploitation of chemical profiles by conjugate variable analysis
M. Gay et al.

Title Page
Abstract Introduction

Conclusions References
Tables Figures

Back Close
Full salars present in Northern Chile and Bolivia, mostly during the dry season.In tropical South America, the fire season takes place during the dry season, and fire events recorded along the Illimani ice core were identified by simultaneous peaks of ammonium, formate and nitrate and significant amounts of soluble potassium (de Angelis et al., 2008).Local emissions of nitrate and ammonium by soils and formate by vegetation may also be expected more likely during the warm wet season.Oxalate is directly emitted through biomass burning, but may also be produced through the oxidation of volatile organic compounds (VOC) scarcely investigated in the South American atmosphere.

Optimization
Optimal values of the high pass filter cut-off frequency f b and the median filter length M are determined for the seven new ions using the criterion presented in Sect.4.2.4 with reference to Pinatubo and Mout Agung eruption dates.The results of the optimization procedure are summarized in Table 2, along with the corresponding minimum value of F (f b , M) defined in Eq. ( 17).We also specify if the 1 yr −1 frequency was found in the corresponding concentration spectra expressed as a function of the time delay estimate (Fig. 14).The eight ions are sorted according to ascending values of the distance criterion F (f b , L).This makes it possible to experimentally estimate the relative quality of chronologies based on each concentration depth profile.Peaks corresponding to annual frequency are clearly visible in nitrate, oxalate, calcium, magnesium and potassium panels in Fig. 14, where concentration spectra are expressed as a function of estimate time for the eight ions.This is not true for sulfate, sodium and formate and is consistent with the sorting presented in Table 2 according

Chronology overall comparison
We are aware that the in-depth understanding of the spectral characteristics of signals extracted from ion concentration depth profiles requires further investigations to answer many questions, which is not beyond the scope of this paper.However, a very preliminary interpretation of signal analysis in terms of time marker reliability may be extracted from Fig. 15, where time-depth relationships respectively deduced from NO − 3 , Ox, Ca 2+ , Mg 2+ and K + concentration profiles are shown along with the time-depth relationship used by de Angelis et al. (2003).A rather good agreement between the five chronologies proposed here as well as with the chronology used by de de Angelis et al. (2003) is observed down to ca. 28 m w.e., i.e. from the late 1950s to the date corresponding to the signal origin depth.Differences increase with depth, reaching a few years at 36 m w.e. for K + and Mg 2+ but with opposite signs.A more detailed representation of chronology trends may be found in Fig. 16 which shows the deviation of our eight independent chronologies from the chronology of de Angelis et al. (2003).Nitrate and calcium provide the best agreement with the chronology published by de Angelis et al. (2003) and the two volcanic tie points.This confirms our hypothesis that these two ions are produced during the dry season: nitrate production is largely dominated by biomass burning compared with soil emission during the warm season and windblown soil dust is by far the most important source of soluble calcium, the relative contribution of volcanic ash remaining low.The good performance of oxalate, both in terms of seasonality and dating accuracy, is very likely due to enhanced production through biomass burning combined with heterogeneous uptake on mineral dust during the dry season as supported by trend similarities in calcium and oxalate deviation curves.Several mechanisms may contribute to overlap the signal seasonality related to emission and/or transport processes of the other ions.Like calcium, magnesium and the part of potassium not related to biomass burning come mainly from erodible soil dust emission, with a sporadic volcanic contribution.However, the composition and dissolution rate of soil dust and volcanic ashes significantly vary depending on soil Introduction

Conclusions References
Tables Figures

Back Close
Full and eruption types influencing the signal of the less soluble species (i.e.magnesium, potassium) in a more significant way than the signal of the more soluble ones like calcium.Episodic transport of polluted plumes of local origin through convection during the wet season since the mid 1950s in addition to the timing of stratospheric volcanic arrivals are very likely responsible for the absence of clear 1 yr −1 frequency in sulfate profile and corresponding ice age overestimate.In a similar way, the occurrence of non periodic peaks related to halite formation and to partial dissolution of ashes in volcanic plumes may become significant compared to seasonal halite emission of sodium and thus explain the bad performance of this ion.Despite a similar origin in relation with biomass emissions the formate signal does not show the 1 yr −1 frequency clearly visible on the nitrate spectrum: formate may be deposited either as salt and then co-vary with nitrogen species, or as formic acid which is known to be sensitive to post deposition processes including peak removal from strongly acidic layers and relocation on mineral particles.

Conclusions
Spectral analysis is considered in this study to be a potential ice core dating tool and is applied for the first time to chemical concentration profiles assumed to be at least partly modulated by seasonal variability of atmospheric content and precipitation patterns.This chronology approach does not require the preliminary knowledge of net accumulation rate or ice flow characteristics and was checked on a tropical Andean ice  (2003).Several results may be highlighted: -Both methods gave the spatial frequency as a function of depth.
-Spectrum replication was detected on both the scalogram and the spectrogram below 48 m w.e. and we demonstrate that replication occured when spatial frequency became greater than the experimental Shannon frequency imposed by the initial chemical sampling procedure.We propose a method of counter balancing aliasing that notably improves ice core chronology in comparison with manual counting of pseudo-annual layer as published by Ramirez et al. (2003).This observation may be useful for glaciologists, allowing them to adapt the sampling resolution at sites where an estimate, even very rough, of annual snow accumulation rate and ice thinning is available.However, it has to be underlined that spectrum replication significantly increased the difficulties encountered here, preventing us from directly assessing the chronology accuracy along the aliased part of the data.
-The chronology provided by the Fourier conjugate analysis for the part of the ice core located above the aliasing depth, i.e. above 48 m w.e., was in close agreement with the chronology based on the multi-proxy year-by-year counting of wet season minima and the Pinatubo and Mont Agung eruptions dates.Wavelets transforms appeared to be less efficient: the scalogram is disrupted by frequent discontinuities related to transients in chemical signals, the seasonal variability of the nitrate signal was not detected, and large discrepancies were observed with the chronology published by de Angelis et al. (2003).
The Fourier conjugate analysis was then applied to the concentration profiles of seven other ions considered as valuable although intricate markers of various biogeochemical reservoirs.Comparing the eight independent chronologies available provided useful information on the natural variability of the corresponding ions.Introduction

Conclusions References
Tables Figures

Back Close
Full Much faster than multi-parameter approaches based on the visual identification of seasonal extrema and not influenced by subjective factors, this innovative chronology approach appears to be particularly suitable for ice cores extracted at poorly documented sites, where the dynamics of seasonal signals may be significantly influenced by source proximity or regional meteorological features, and where emissions of local or regional origin make it difficult to unambiguously identify particular events of global importance or anthropogenic trends generally used as adjusting marks along polar ice cores.When applied to chemical records extracted from polar deep ice cores, the Fourier conjugate analysis appears as a useful tool for studying the variability of complex atmospheric phenomena at annual as well as longer time scale.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |spatial frequency is corroborated by the two autocorrelation functions (seeWheatley et al., 2012)  shown in Fig.2.These autocorrelation functions are calculated on 128 samples.The red circles depict the autocorrelation function calculated between 14.92 m w.e. and 19.7 m w.e., and the blue dots the autocorrelation function calculated between 22.4 m w.e. and 27.15 m w.e.The first maximum of every autocorrelation function provides an estimate of annual layer thickness, i.e.: 0.59 m w.e. as a mean between 14.92 m w.e. and 19.7 m w.e.(red circles), that implies a spatial frequency of 1.69 m w.e.−1 and 0.41 m w.e. as a mean between 22.4 m w.e. and 27.15 m w.e.(blue dots), that implies a spatial frequency of 2.46 m w.e.−1 .Annual layer thinning with depth is related to firn densification and ice flowing.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | We have a set of 1697 samples of c NO − 3 (n).The nitrate concentration function c NO − 3 Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | where h(d ) is a sliding window of length N (Hanning function).The layer thickness resolution of the spectrogram is determined by the thickness of the selected sliding window h(d ).
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Spatial frequency estimates from the spectrogram and the scalogram are given by the spatial frequency f s of the spectrogram and scalogram maxima at each depth, d : Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 9
Figure 9 shows the spatial frequency estimated with the spectrogram (left) and the scalogram (right) as a function of depth and obtained with the optimum values of f b and M. Discontinuities are smoothed on the two plots.Their attenuation is more efficient on the spectrogram estimate than on the scalogram estimate of the spatial frequency.
8 and 9 comes from aliasing.The aliasing is corrected by replacing the decreasing part of the spectrum by its symetrical representation with respect to the maximum frequency observed just before aliasing.Maximum frequency values observed before aliasing are f s,spec,max = 8.59 (m w.e.) −1 and f s,scal,max = 8.37 (m w.e.) −1 for the spectrogram and the scalogram estimates, respectively.Then for depths d greater than the aliasing depth d aliasing , we replace f s (d ) by (2f smax − f s (d )).The resulting estimates of spatial frequencies are presented in Fig. 10 (left panel: spectrogram; right panel: scalogram).The trend of spatial frequency estimates increases with depth, which is coherent with firn densification.corresponding to the depth d as a function of the spatial frequency f s (u), of the signal analysis depth origin d min , and of the date t(d min ) associated with d min .The first measurement based on spectrogram window of length N and a median filter window of length M is related to the depth ∆(d ): 128, M = 65 and P s = 0.0373 m w.e., the depth of this first measurement d min is: d min = 3.6 m we (22) We estimate t(d min) by t(d min ), assuming that f s (u) is constant on the segment [0, d min ]: t(d min ) = f (d min ) • d min (23) This leads to t(d min ) = 5.28 yr for the spectrogram and t(d min ) = 5.33 yr for the scalogram.Finally, the estimates, t(d ), of the time delay between the drilling date and the CPD Discussion Paper | Discussion Paper | Discussion Paper | and the time estimates t(d ) given by the two chronologies respectively resulting from Fourier conjugate analysis and wavelet transforms and using a spline transformation we get the concentration of NO − 3 as a function of the time estimate: S NO − Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | time scale published in 2003 for the upper part of the core (dating 1) is represented by a black line in the upper part of Fig. 13 (top panels) while time scales respectively deduced from the S NO − 3 spectrogram (left) and scalogram (right) are represented by a green line (dating 2).The corresponding differences are presented in the bottom part of the figure.The date gap corresponding to the spectrogram (left panel) remains close to zero from the signal analysis origin depth to 30 m w.e.It slightly increases deeper down, reaching 1.5 yr at 35 m w.e.The chronology deviation corresponding to the scalogram (right panel) is higher, particularly between 5 and 15 m w.e. and reaches 4 yr at 10 m w.e.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | core recovered in the summit area of the Nevado Illimani (6350 m a.s.l.).Fourier conjugate analysis and orthogonal Morlet wavelet transforms, traditionally used to study long term climate variability once core chronology is achieved, were first compared for tracking spatial frequency in nitrate concentration depth profile.Fourier conjugate analysis uses a spectrogram, while Morlet wavelet transforms use a scalogram.High pass and median filters aimed at correcting spectrogram and scalogram discontinuities were optimized using two volcanic layers attributed to the Pinatubo (June 1991) and Discussion Paper | Discussion Paper | Discussion Paper | Mount Agung (March 1963) eruptions by de Angelis et al.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Table 1 .
Dating of tie points from the spectrogram and the scalogram.