Milankovic Pseudocycles Recorded in Sediments and Ice Cores Extracted by Singular Spectrum Analysis
 ^{1}Institut de Physique du Globe, Paris University, Paris 75005, France
 ^{2}Sorbonne Université, CNRS, EPHE, UMR 7619 METIS, 75005 Paris, France
 ^{3}Genome Structure and Instability Laboratory, CNRS UMR 7196, INSERM U1154, Muséum national d’Histoire naturelle, Alliance Sorbonne Université, 75005 Paris, France
 ^{4}Centre de Microscopie et d’Imagerie Numérique, Muséum National d’Histoire Naturelle, Paris, France, Muséum national d’Histoire naturelle, Alliance Sorbonne Université, 75005 Paris, France
 ^{1}Institut de Physique du Globe, Paris University, Paris 75005, France
 ^{2}Sorbonne Université, CNRS, EPHE, UMR 7619 METIS, 75005 Paris, France
 ^{3}Genome Structure and Instability Laboratory, CNRS UMR 7196, INSERM U1154, Muséum national d’Histoire naturelle, Alliance Sorbonne Université, 75005 Paris, France
 ^{4}Centre de Microscopie et d’Imagerie Numérique, Muséum National d’Histoire Naturelle, Paris, France, Muséum national d’Histoire naturelle, Alliance Sorbonne Université, 75005 Paris, France
Abstract. Milankovic cycles describe the changes in the Earth's orbit and rotation axis and their impact on its climate over thousands of years. Singular Spectrum Analysis (SSA) is a signal processing method that is best known for its ability to find and extract pseudocycles in complex signals. In this short paper, we propose to apply it to three time series that have been proposed as geological reference time scales, in order to retrieve, compare and identify their Milankovic periodicities: (1) LR04, a stack of PlioPleistocene benthic microfossil records (Lisiecki and Raymo, 2005), (2) the CO_{2} and CH_{4} records from the Vostok ice core (Petit et al, 1999) and (3) the longterm orbital solution La04 for the insolation of Laskar et al (2004). The Vostok CO_{2} and CH_{4} series share the first 7 SSA components, three main ones (98, 104, 39 kyr), and four smaller ones (18, 22, 65, 180 kyr). CO_{2} displays a component at 28 kyr and a doublet at 61 and 62 kyr. CH_{4} displays a doublet near 50 kyr. 18/22 ky is a precession doublet, 62 kyr an insolation component, and 95/105 kyr an insolation/eccentricity doublet. The 49/50 kyr doublet in CH_{4} is not found in the orbital model. The SSA results for the La04 orbital solution are in excellent agreement with the values obtained by Laskar et al (2004). Four SSA components of obliquity are almost identical (rounded figures are 41, 54, 29 and 39 kyr). As far as eccentricity is concerned, the first five components are 404, 95, 124, 99, and 132 kyr. The next components are not found in our list of components for eccentricity, but they are in the SSA of insolation, at 2338, 970, 488 and 684 kyr. With more than 20 components, the LR04 stack is the richest series. In order of decreasing amplitude, one encounters 41, 95 and 75 kyr components. Next are smaller 39.5 and 53.6 kyr components, and a 22.4 kyr component. One recognizes one of the two main precession components, the doublet of obliquity components, a line at 47.4 kyr that is not found in any of the other spectra, and a doublet at 53.6 and 55.7 kyr, corresponding to the line at 54 kyr found in all four orbital quantities. Next comes a line at 63.6 kyr that may correspond to a line in insolation, CH_{4} and CO_{2}. Then come components from eccentricity variations at 75.2, 94.5, 107.2, 132.1, 198.6 and 400.9 kyr. The remaining components of LR04 show up in La04. The “elusive ~200 kyr eccentricity cycle” of Hilgen et al (2020) is actually present in all three series, in the La04 orbital model as a 195 ± 6 kyr component of eccentricity and in LR04 as a 198.6 ± 5.6 kyr component. Finding not only the main expected Milankovic periodicities but also many “secondary” components with much smaller amplitudes gives confidence in our iterative SSA method (iSSA), on the quality of the La04 model and on the remarkable LR04 sedimentary stack, with more than 15 “ Milankovic periods”.
Fernando Lopes et al.
Status: closed

RC1: 'Milankovitch forcing is about SUMMER insolation, not global !  Detection is about significance, not profusion !', Anonymous Referee #1, 20 Dec 2021
I have difficulties to understand the main topic of this manuscript: Either the main message is that paleoclimatic records contain astronomical periodicities, something which is wellknown for almost 50 years. Or, more probably, the aim is to present a “new” spectral analysis method that could detect extremely faint periodicities in geological records. In the first case, I would reject the paper on the ground that it brings nothing new. But the authors are not even citing the famous foundational Hays et al (1976) paper in their reference list, which is, to say the least, very peculiar and which suggest that they may not be familiar with this topic and its extremely abundant literature. Besides (see below), it appears that the authors have even missed the main point of Milankovitch theory, which is based on summer insolation, and not global mean insolation... In the second case, I would also reject the paper since it does not show any information on the key question of spectral analysis: statistical relevance. Indeed, the authors seem to be unaware that white noise contains “all frequencies”. Therefore, finding frequencies in a signal is NOT (has never been) the key point of spectral analysis. The only interesting question is to estimate their statistical significance. But there is simply no mention at all of statistics in this paper. On these two main grounds, and several others listed below, I do not recommend publication of this manuscript.
Major comments on SSA.
SSA is now a very classical spectral analysis method. It was largely developed in the 1980s (Broomhead and King, 1986) and was used for paleoclimatic studies since then (Vautard and Ghil, 1989). It was even applied to the Vostok record as early as 1994 (Yiou et al., 1994). Again, it is strange not to find these references in a scientific paper on SSA applied to paleoclimate records in general and to the Vostok data in particular. In fact, SSA has been applied to Vostok and to LR04 in many previous papers over the last 20 years : it is so classical that this analysis is now even part of some student textbooks (eg. Martinson 2018, page 531, Figure 15.6 showing the SSA analysis of LR04). In any case, it is difficult to find any novelty in this paper.
Spectral analysis is a branch of statistics. Finding many frequencies (real or spurious) in a signal is known to be extremely easy, since the discovery of the Fourier transform. In particular, finding astronomical frequencies and their numerous harmonics in paleoclimatic signals is certainly not surprising. Actually, the opposite would be much more unexpected. This is even more true for the insolation itself !… Not finding the main and secondary astronomical periodicities that are used as input of the insolation computation would be a sure signature of some methodological error. I therefore do not understand the value of applying a spectral analysis method (in particular on insolation), without discussing the "detection level" problem. Indeed, the key mathematical difficulty is to decide whether it is “noise” or “signal” using some statistical test. There is no mention of any statistics in this paper, which certainly disqualifies it entirely. This is strange since standard SSA tools used in the climate community (eg. SSAMTM toolkit) do provide some statistical tools and safeguards. The harmonics and combination tones found by the authors have often a very small amplitude (below 1% or 0.1% of the variance) and to prove their significance is certainly not an easy task, given that the records are rather “short” and “noisy”. This is precisely the main difficulty of spectral analysis, and a problem that climate scientists are usually facing. But the authors of this paper are not even mentioning or discussing this central point, as if they were unaware of the main problem of “spectral analysis estimation”.
In addition, it is worth mentioning that detecting several significant frequencies (let alone NUMEROUS significant frequencies) is a much more difficult problem than detecting a single one. Indeed, the nullhypothesis cannot be a simple (red or white) noise hypothesis. A classical workaround is to perform MonteCarlo simulations of synthetic signals with and without each single periodicity and assess whether the results are statistically distinguishable. But this would probably be a very different paper.
Just to be more precise: Of course, I do not mention the uncertainty in the periodicities as listed in table 1 (obtained after SSA decomposition and resampling of the data), but the significance of the SSA components themselves. In particular, the KEY parameter here is the choice of the embedding dimension. If chosen too high (extremely likely here), then there are a lot of SSA components with dubious relevance. The same phenomenon applies to all spectral methods: the more degrees of freedom you are allowing, the more frequencies you get, but the less significant they are. What is the embedding dimension of the analysis here? This is not even mentionned ! But according to line 62, “SSA was able to extract up to 75 cycles”. My guess is that the embedding dimension is a few hundreds, which represents quite a lot of freedom for the method to generate a lot of “spurious” frequencies.
It is very strange to use the old Vostok record (Petit et al. 1999) which covers only 4 glacial cycles while a more recent one (Dome C), covering twice as many cycles, has been available for more than a decade now (see eg: Guo et al (2012) for the SSA analysis of CH4 at Dome C over 8 cycles, not 4 as here). Obviously, the results are likely to be statistically more robust with a longer time series, if statistics matters of course.
“Finding not only the main expected Milankovic periodicities but also many “secondary” components with much smaller amplitudes gives confidence in our iterative SSA method”.
But I would expect standard methods to do the same, with attributing to these periodicities a very low confidence. In any case again, finding many frequencies is easy and NOT the main problem of “spectral analysis”. Besides, what is the subject of the paper? Is it about climate or mathematics?
Major comments on Climate
The authors appear unaware that the word “insolation” is not sufficient to describe the astronomical forcing of climate. They are presenting the spectral analysis of “the insolation of Laskar et al (2004)”. This is clearly a very insufficient description. When looking more closely at the spectrum (Figure 3), it seems that they are in fact discussing global averaged insolation (linked to eccentricity only)! Otherwise the spectrum would be dominated by the 41k and the 23k periodicities. Since there is no power in the obliquity and precession band, the only explanation is that the authors are using global averaged insolation. But this is a major misunderstanding on the Milankovitch forcing. Indeed, it is well known since the 19^{th} century that global insolation changes are much too weak to affect climate and scientists have discussed local seasonal insolation, not mean global one. When climate scientists talk about “insolation” in the context of glacial cycles, they are talking about summer insolation at high northern latitude according to Milankovitch theory. THIS IS NOT what is used in this analysis. In other words, this paper has missed entirely its claimed goal of discussing the Milankovitch forcing, simply by looking at the wrong insolation…
To conclude, its seems to me that the authors have misunderstood both the aim of "spectral analysis" and the foundations of the "Milankovitch theory". I cannot be favorable to publication.
Some references
Broomhead D.S., King G.P. (1986) Extracting qualitative dynamics from experimental data. Physica D 20: 217236.
Vautard R., Ghil M (1989) Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. Physica D 35: 295424.
Yiou, P., Ghil, M., Jouzel, J., Paillard, D., & Vautard, R. (1994). Nonlinear variability of the climatic system, from singular and power spectra of Late Quaternary Records. Climate Dynamics, 9, 371–389.
Guo Z, Zhou, X., Wu H. (2012) Glacialinterglacial water cycle, global monsoon and atmospheric methane changes, Clim. Dyn 39, 10731092.
Viaggi P. (2018) δ^{18}O and SST signal decomposition and dynamic of the PliocenePleistocene climate system: new insights on orbital nonlinear behavior vs. longterm trend, Progress in Earth an Planetary Science.
Martinson, D. (2018). Empirical Orthogonal Function (EOF) Analysis. In Quantitative Methods of Data Analysis for the Physical Sciences and Engineering (pp. 495534). Cambridge: Cambridge University Press. doi:10.1017/9781139342568.016

CC1: 'Reply on RC1', Fernando Lopes, 01 Feb 2022
The comment was uploaded in the form of a supplement: https://cp.copernicus.org/preprints/cp2021126/cp2021126CC1supplement.pdf

CC1: 'Reply on RC1', Fernando Lopes, 01 Feb 2022

RC2: 'Comment on cp2021126', Anonymous Referee #2, 09 Feb 2022
Initial note: I should emphasize that I read the comments of reviewer #1, as well as the authors' rebuttal, only after drafting my own report (see below). From both documents, I am strongly concerned whether the authors' own opinion on the scope of their analysis and the type of information that this analysis can provide to the paleoclimate community is compatible with that of the vast majority of climate scientists and, thus, the readership of Climate of the Past.
Despite I do not necessarily agree with all details of the report of reviewer #1, I share their opinion that SSA is used here as a statistical approach, even if the authors strictly oppose to this view. (In the end, the decomposition is related to autocorrelation coefficients unfolded into a matrix, i.e., a basic statistical quantity that is estimated from the data under study.) Hence, I agree and insist that the findings reported need to be evaluated using standard scientific measures applying to any kind of statistical method, i.e., combining a descriptive statistics (identified oscillation modes) with an associated significance assessment. The authors argue in their rebuttal against this view in a way that is worrying me to a large extent, since they appear to miss the essential points of criticism raised by our colleague.
This being said, let me come to the results of my own reading of this manuscipt:
_____
The authors present an analysis of periodic components embedded in several paleoclimate records, namely different proxies from the Antarctic Vostok ice core record and the LR04 PlioPleistocene benthic oxygen isotope stack, together with some insolation variations and orbital cycles during the last 250 Myr. For this purpose, they utilize a variant of the celebrated singular spectrum analysis to decompose the different time series based on lagged trajectory matrices and associate the corresponding pairs of eigenvectors(?) with oscillatory modes with welldefined frequency content. The results are then compared among the different time series to highlight the presence of certain orbital cycles and combinations thereof in the paleoclimate records.
From their publication history, the authors appear to be well familiar with the used statistical methodology, which is probably more than they should expect from the average readership of Climate of the Past. They may have the feeling that it is sufficient to introduce the iSSA technique by referring to their expressive own publication list; I personally do not think that this is appropriate. To familiarize the readership with the specifications of the method, the method needs to be thoroughly explained as part of the manuscript. I may belong to a minority knowing well about SSA in general and even having it applied in my own previous work, but I am confident that the vast majority of possible readers will not have this background. However, even I do not know what exactly makes up the special feature of the used Iterative SSA without further explanation. The authors are kindly requested to add more details on the statistical methodology to familiarize their readers with the used approach.
Somewhat related to this, my understanding of SSA has always been that it allows to identify periodic components in time series without the normal restriction of Fourier analysis assuming a harmonic shape of oscillations. However, to my best knowledge, it still requires a fixed repetitive pattern (especially with fixed oscillation frequency) for each oscillatory mode. If this is correct, I am worrying about justifying the known nonstationarity of Milankovich components in paleoclimate records, which is well visible in the LR04 dataset in terms of ice age cycles changing from approximately 41 kyr to 100 kyr along with the midPleistocene transition, as barely more than amplitude changes of those (nonharmonic) oscillations. A more explicit discussion of this aspect (along with the restriction to oscillations with a fixed frequency) appears essential for making the presented analysis useful for the paleoclimate community. In essence, this point is important for the motivation of using (i)SSA in this study. Why not employing any other time series decomposition technique that might more flexibly cope with nonstationarity not only in amplitude, but also frequency? The authors motivate the use of iSSA exclusively by the fact that they have used this method successfully to other datasets in the past – in my opinion, a proper motivation would need to start with some (justifiable) working hypothesis on the type of oscillations to be identified.
Even when neglecting the aforementioned issues, I would still need some clarification about the identification of the different oscillatory modes. I suppose that they correspond to pairs of SSA eigenvectors/singular vectors that are in quadrature, yet many of them do likely not only lack physical origin, but may also have low magnitude making them practically indistinguishable from noise. From a statistical perspective, I may argue that one should only be interested in components whose variability amplitude clearly exceeds a certain significance threshold. However, I did not find a word about any significance assessment related to the reported oscillatory modes. If there has been some corresponding assessment that was simply not (yet) reported: what has been the underlying null model? (White noise, correlated noise, etc.?)
Finally, if I accept that the list of oscillatory modes identified for each time series is useful for some scientific purpose and considers certain significance criteria, I still wonder about the interpretation. In this regard, I found the wording/terminology used by the authors pretty vague or even somewhat confusing and inconsistent – periodicity, cycle, pseudocycle, quasiperiodicity, etc. Some of these words appear to have a precise meaning (and especially “quasiperiodic” typically means something distinctively different than what it is used for here), others are more colloquial. I would warmly welcome a more “mechanistic” discussion on which of the cycles are actually “relevant” and of physical origin, which rather occur due to (nonlinear?) combinations of other (nonharmonic) “cycles”, and which may be barely more than statistical artifacts (e.g., affected by the time step of the considered series). Uncovering and explaining this thoroughly would be a strong contribution that could actually be helpful for the readers of this manuscript and the whole community. To this end, I found every “attempt” of providing corresponding hints in this work very vague, to say the least.
In summary, I am not convinced that the paper as it stands now should be published in Climate of the Past. I have strong sympathy for this type of analysis, but the way it is reported here it is very technical, not very appropriate for the readership of the journal, and leaves open many important aspects.
Specific comments:
The abstract should reflect the main messages, not a detailed enumeration of all specific results. Please focus on the main insights of this study.
Line 15: What do you define as a “Milankovic periodicity”? Typically, I would associate this term only with the classical major oscillation frequencies of the Earth’s orbital parameters, while the authors seem to generously expand this term to everything they find in their analysis.
The Vostok and LR04 data have undergone a sequence of more or less heavy preprocessing steps (interpolation, stacking) that may have modified their spectral characteristics especially in the higherfrequency part beyond simple age uncertainty effects (e.g., any interpolation to a different time axis necessarily induces correlations). This important source of uncertainty (and potentially error for the presented analysis) needs to be discussed.
The authors report that LR04 is the richest series in terms of the number of embedded oscillatory components. To me, this is not quite surprising: other than the orbital solution by Laskar, the series is not strictly deterministic and not smooth, which likely introduces additional components into the decomposition “just by chance”. Moreover, since it is considerably longer than the Vostok core, it can also contain more lowerfrequency modes than the latter. Both features are in my opinion well traceable in Tab. 1.
What is the insolation time series used by the authors? Summer insolation at 65°N? Please specify!
Line 43: What does plate tectonics have to do with magnetic reversals?
Line 80: What do you mean by “ratio planktonic to ice”?
Figure 2: Why do you use a reverted time axis as compared to Fig. 1?
Line 104: “but not for atmospheric” – atmospheric what?
Line 119: Does Laskar et al. 2004 really use nine planets? Rather eight planets? (Plus Pluto? Other large bodies with larger gravitational effect on the Earth’s orbit than this one?)
Lines 135144: I do not see any relevance of this paragraph for the present analysis, despite elevating the number of selfcitations.
Generally, the wording of “a line” for referring to the frequency of any SSA mode is not quite useful.
Line 281282: This sentence provides an important hint on the nature of the secondary/pseudocycles found in this analysis. This should be further detailed.
Line 302: The phase of an oscillation is continuously varying, so this does not need to be mentioned. The fact that the amplitudes of SSA modes can vary with time is important for the present analysis, but the corresponding variations are not really exploited in detail for the reported analysis. What seems to be relevant to me are slow modulations of the amplitudes, which could indicate some kind of crossfrequency coupling between different modes/processes.
Line 308: What do you mean by “response functions of the various components”?
Lines 336343: This is not a result of the present work and therefore completely irrelevant for the conclusions.
In general, the conclusions section leaves the same impression as the abstract, being full of technical details making the reader lost with trying to identify the key findings of the paper. In any case, this part reads very repetitive.
Can the differences between the Fourier analysis of Laskar et al. and the iSSA results for the same time series reported in this work be simply explained by the fact that one method explicitly assumes harmonic oscillations while the other does not?
Line 371: “Eccentricity has the main influence on insolation.” I would strongly doubt that the physical effect (in W/m2) of eccentricity does actually exceed those of the other orbital components.
Lines 384386: I think it is statistically meaningless to discuss apparent cyclic components with long periods that only fit two or three times into the total length of the studied time series.
Status: closed

RC1: 'Milankovitch forcing is about SUMMER insolation, not global !  Detection is about significance, not profusion !', Anonymous Referee #1, 20 Dec 2021
I have difficulties to understand the main topic of this manuscript: Either the main message is that paleoclimatic records contain astronomical periodicities, something which is wellknown for almost 50 years. Or, more probably, the aim is to present a “new” spectral analysis method that could detect extremely faint periodicities in geological records. In the first case, I would reject the paper on the ground that it brings nothing new. But the authors are not even citing the famous foundational Hays et al (1976) paper in their reference list, which is, to say the least, very peculiar and which suggest that they may not be familiar with this topic and its extremely abundant literature. Besides (see below), it appears that the authors have even missed the main point of Milankovitch theory, which is based on summer insolation, and not global mean insolation... In the second case, I would also reject the paper since it does not show any information on the key question of spectral analysis: statistical relevance. Indeed, the authors seem to be unaware that white noise contains “all frequencies”. Therefore, finding frequencies in a signal is NOT (has never been) the key point of spectral analysis. The only interesting question is to estimate their statistical significance. But there is simply no mention at all of statistics in this paper. On these two main grounds, and several others listed below, I do not recommend publication of this manuscript.
Major comments on SSA.
SSA is now a very classical spectral analysis method. It was largely developed in the 1980s (Broomhead and King, 1986) and was used for paleoclimatic studies since then (Vautard and Ghil, 1989). It was even applied to the Vostok record as early as 1994 (Yiou et al., 1994). Again, it is strange not to find these references in a scientific paper on SSA applied to paleoclimate records in general and to the Vostok data in particular. In fact, SSA has been applied to Vostok and to LR04 in many previous papers over the last 20 years : it is so classical that this analysis is now even part of some student textbooks (eg. Martinson 2018, page 531, Figure 15.6 showing the SSA analysis of LR04). In any case, it is difficult to find any novelty in this paper.
Spectral analysis is a branch of statistics. Finding many frequencies (real or spurious) in a signal is known to be extremely easy, since the discovery of the Fourier transform. In particular, finding astronomical frequencies and their numerous harmonics in paleoclimatic signals is certainly not surprising. Actually, the opposite would be much more unexpected. This is even more true for the insolation itself !… Not finding the main and secondary astronomical periodicities that are used as input of the insolation computation would be a sure signature of some methodological error. I therefore do not understand the value of applying a spectral analysis method (in particular on insolation), without discussing the "detection level" problem. Indeed, the key mathematical difficulty is to decide whether it is “noise” or “signal” using some statistical test. There is no mention of any statistics in this paper, which certainly disqualifies it entirely. This is strange since standard SSA tools used in the climate community (eg. SSAMTM toolkit) do provide some statistical tools and safeguards. The harmonics and combination tones found by the authors have often a very small amplitude (below 1% or 0.1% of the variance) and to prove their significance is certainly not an easy task, given that the records are rather “short” and “noisy”. This is precisely the main difficulty of spectral analysis, and a problem that climate scientists are usually facing. But the authors of this paper are not even mentioning or discussing this central point, as if they were unaware of the main problem of “spectral analysis estimation”.
In addition, it is worth mentioning that detecting several significant frequencies (let alone NUMEROUS significant frequencies) is a much more difficult problem than detecting a single one. Indeed, the nullhypothesis cannot be a simple (red or white) noise hypothesis. A classical workaround is to perform MonteCarlo simulations of synthetic signals with and without each single periodicity and assess whether the results are statistically distinguishable. But this would probably be a very different paper.
Just to be more precise: Of course, I do not mention the uncertainty in the periodicities as listed in table 1 (obtained after SSA decomposition and resampling of the data), but the significance of the SSA components themselves. In particular, the KEY parameter here is the choice of the embedding dimension. If chosen too high (extremely likely here), then there are a lot of SSA components with dubious relevance. The same phenomenon applies to all spectral methods: the more degrees of freedom you are allowing, the more frequencies you get, but the less significant they are. What is the embedding dimension of the analysis here? This is not even mentionned ! But according to line 62, “SSA was able to extract up to 75 cycles”. My guess is that the embedding dimension is a few hundreds, which represents quite a lot of freedom for the method to generate a lot of “spurious” frequencies.
It is very strange to use the old Vostok record (Petit et al. 1999) which covers only 4 glacial cycles while a more recent one (Dome C), covering twice as many cycles, has been available for more than a decade now (see eg: Guo et al (2012) for the SSA analysis of CH4 at Dome C over 8 cycles, not 4 as here). Obviously, the results are likely to be statistically more robust with a longer time series, if statistics matters of course.
“Finding not only the main expected Milankovic periodicities but also many “secondary” components with much smaller amplitudes gives confidence in our iterative SSA method”.
But I would expect standard methods to do the same, with attributing to these periodicities a very low confidence. In any case again, finding many frequencies is easy and NOT the main problem of “spectral analysis”. Besides, what is the subject of the paper? Is it about climate or mathematics?
Major comments on Climate
The authors appear unaware that the word “insolation” is not sufficient to describe the astronomical forcing of climate. They are presenting the spectral analysis of “the insolation of Laskar et al (2004)”. This is clearly a very insufficient description. When looking more closely at the spectrum (Figure 3), it seems that they are in fact discussing global averaged insolation (linked to eccentricity only)! Otherwise the spectrum would be dominated by the 41k and the 23k periodicities. Since there is no power in the obliquity and precession band, the only explanation is that the authors are using global averaged insolation. But this is a major misunderstanding on the Milankovitch forcing. Indeed, it is well known since the 19^{th} century that global insolation changes are much too weak to affect climate and scientists have discussed local seasonal insolation, not mean global one. When climate scientists talk about “insolation” in the context of glacial cycles, they are talking about summer insolation at high northern latitude according to Milankovitch theory. THIS IS NOT what is used in this analysis. In other words, this paper has missed entirely its claimed goal of discussing the Milankovitch forcing, simply by looking at the wrong insolation…
To conclude, its seems to me that the authors have misunderstood both the aim of "spectral analysis" and the foundations of the "Milankovitch theory". I cannot be favorable to publication.
Some references
Broomhead D.S., King G.P. (1986) Extracting qualitative dynamics from experimental data. Physica D 20: 217236.
Vautard R., Ghil M (1989) Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. Physica D 35: 295424.
Yiou, P., Ghil, M., Jouzel, J., Paillard, D., & Vautard, R. (1994). Nonlinear variability of the climatic system, from singular and power spectra of Late Quaternary Records. Climate Dynamics, 9, 371–389.
Guo Z, Zhou, X., Wu H. (2012) Glacialinterglacial water cycle, global monsoon and atmospheric methane changes, Clim. Dyn 39, 10731092.
Viaggi P. (2018) δ^{18}O and SST signal decomposition and dynamic of the PliocenePleistocene climate system: new insights on orbital nonlinear behavior vs. longterm trend, Progress in Earth an Planetary Science.
Martinson, D. (2018). Empirical Orthogonal Function (EOF) Analysis. In Quantitative Methods of Data Analysis for the Physical Sciences and Engineering (pp. 495534). Cambridge: Cambridge University Press. doi:10.1017/9781139342568.016

CC1: 'Reply on RC1', Fernando Lopes, 01 Feb 2022
The comment was uploaded in the form of a supplement: https://cp.copernicus.org/preprints/cp2021126/cp2021126CC1supplement.pdf

CC1: 'Reply on RC1', Fernando Lopes, 01 Feb 2022

RC2: 'Comment on cp2021126', Anonymous Referee #2, 09 Feb 2022
Initial note: I should emphasize that I read the comments of reviewer #1, as well as the authors' rebuttal, only after drafting my own report (see below). From both documents, I am strongly concerned whether the authors' own opinion on the scope of their analysis and the type of information that this analysis can provide to the paleoclimate community is compatible with that of the vast majority of climate scientists and, thus, the readership of Climate of the Past.
Despite I do not necessarily agree with all details of the report of reviewer #1, I share their opinion that SSA is used here as a statistical approach, even if the authors strictly oppose to this view. (In the end, the decomposition is related to autocorrelation coefficients unfolded into a matrix, i.e., a basic statistical quantity that is estimated from the data under study.) Hence, I agree and insist that the findings reported need to be evaluated using standard scientific measures applying to any kind of statistical method, i.e., combining a descriptive statistics (identified oscillation modes) with an associated significance assessment. The authors argue in their rebuttal against this view in a way that is worrying me to a large extent, since they appear to miss the essential points of criticism raised by our colleague.
This being said, let me come to the results of my own reading of this manuscipt:
_____
The authors present an analysis of periodic components embedded in several paleoclimate records, namely different proxies from the Antarctic Vostok ice core record and the LR04 PlioPleistocene benthic oxygen isotope stack, together with some insolation variations and orbital cycles during the last 250 Myr. For this purpose, they utilize a variant of the celebrated singular spectrum analysis to decompose the different time series based on lagged trajectory matrices and associate the corresponding pairs of eigenvectors(?) with oscillatory modes with welldefined frequency content. The results are then compared among the different time series to highlight the presence of certain orbital cycles and combinations thereof in the paleoclimate records.
From their publication history, the authors appear to be well familiar with the used statistical methodology, which is probably more than they should expect from the average readership of Climate of the Past. They may have the feeling that it is sufficient to introduce the iSSA technique by referring to their expressive own publication list; I personally do not think that this is appropriate. To familiarize the readership with the specifications of the method, the method needs to be thoroughly explained as part of the manuscript. I may belong to a minority knowing well about SSA in general and even having it applied in my own previous work, but I am confident that the vast majority of possible readers will not have this background. However, even I do not know what exactly makes up the special feature of the used Iterative SSA without further explanation. The authors are kindly requested to add more details on the statistical methodology to familiarize their readers with the used approach.
Somewhat related to this, my understanding of SSA has always been that it allows to identify periodic components in time series without the normal restriction of Fourier analysis assuming a harmonic shape of oscillations. However, to my best knowledge, it still requires a fixed repetitive pattern (especially with fixed oscillation frequency) for each oscillatory mode. If this is correct, I am worrying about justifying the known nonstationarity of Milankovich components in paleoclimate records, which is well visible in the LR04 dataset in terms of ice age cycles changing from approximately 41 kyr to 100 kyr along with the midPleistocene transition, as barely more than amplitude changes of those (nonharmonic) oscillations. A more explicit discussion of this aspect (along with the restriction to oscillations with a fixed frequency) appears essential for making the presented analysis useful for the paleoclimate community. In essence, this point is important for the motivation of using (i)SSA in this study. Why not employing any other time series decomposition technique that might more flexibly cope with nonstationarity not only in amplitude, but also frequency? The authors motivate the use of iSSA exclusively by the fact that they have used this method successfully to other datasets in the past – in my opinion, a proper motivation would need to start with some (justifiable) working hypothesis on the type of oscillations to be identified.
Even when neglecting the aforementioned issues, I would still need some clarification about the identification of the different oscillatory modes. I suppose that they correspond to pairs of SSA eigenvectors/singular vectors that are in quadrature, yet many of them do likely not only lack physical origin, but may also have low magnitude making them practically indistinguishable from noise. From a statistical perspective, I may argue that one should only be interested in components whose variability amplitude clearly exceeds a certain significance threshold. However, I did not find a word about any significance assessment related to the reported oscillatory modes. If there has been some corresponding assessment that was simply not (yet) reported: what has been the underlying null model? (White noise, correlated noise, etc.?)
Finally, if I accept that the list of oscillatory modes identified for each time series is useful for some scientific purpose and considers certain significance criteria, I still wonder about the interpretation. In this regard, I found the wording/terminology used by the authors pretty vague or even somewhat confusing and inconsistent – periodicity, cycle, pseudocycle, quasiperiodicity, etc. Some of these words appear to have a precise meaning (and especially “quasiperiodic” typically means something distinctively different than what it is used for here), others are more colloquial. I would warmly welcome a more “mechanistic” discussion on which of the cycles are actually “relevant” and of physical origin, which rather occur due to (nonlinear?) combinations of other (nonharmonic) “cycles”, and which may be barely more than statistical artifacts (e.g., affected by the time step of the considered series). Uncovering and explaining this thoroughly would be a strong contribution that could actually be helpful for the readers of this manuscript and the whole community. To this end, I found every “attempt” of providing corresponding hints in this work very vague, to say the least.
In summary, I am not convinced that the paper as it stands now should be published in Climate of the Past. I have strong sympathy for this type of analysis, but the way it is reported here it is very technical, not very appropriate for the readership of the journal, and leaves open many important aspects.
Specific comments:
The abstract should reflect the main messages, not a detailed enumeration of all specific results. Please focus on the main insights of this study.
Line 15: What do you define as a “Milankovic periodicity”? Typically, I would associate this term only with the classical major oscillation frequencies of the Earth’s orbital parameters, while the authors seem to generously expand this term to everything they find in their analysis.
The Vostok and LR04 data have undergone a sequence of more or less heavy preprocessing steps (interpolation, stacking) that may have modified their spectral characteristics especially in the higherfrequency part beyond simple age uncertainty effects (e.g., any interpolation to a different time axis necessarily induces correlations). This important source of uncertainty (and potentially error for the presented analysis) needs to be discussed.
The authors report that LR04 is the richest series in terms of the number of embedded oscillatory components. To me, this is not quite surprising: other than the orbital solution by Laskar, the series is not strictly deterministic and not smooth, which likely introduces additional components into the decomposition “just by chance”. Moreover, since it is considerably longer than the Vostok core, it can also contain more lowerfrequency modes than the latter. Both features are in my opinion well traceable in Tab. 1.
What is the insolation time series used by the authors? Summer insolation at 65°N? Please specify!
Line 43: What does plate tectonics have to do with magnetic reversals?
Line 80: What do you mean by “ratio planktonic to ice”?
Figure 2: Why do you use a reverted time axis as compared to Fig. 1?
Line 104: “but not for atmospheric” – atmospheric what?
Line 119: Does Laskar et al. 2004 really use nine planets? Rather eight planets? (Plus Pluto? Other large bodies with larger gravitational effect on the Earth’s orbit than this one?)
Lines 135144: I do not see any relevance of this paragraph for the present analysis, despite elevating the number of selfcitations.
Generally, the wording of “a line” for referring to the frequency of any SSA mode is not quite useful.
Line 281282: This sentence provides an important hint on the nature of the secondary/pseudocycles found in this analysis. This should be further detailed.
Line 302: The phase of an oscillation is continuously varying, so this does not need to be mentioned. The fact that the amplitudes of SSA modes can vary with time is important for the present analysis, but the corresponding variations are not really exploited in detail for the reported analysis. What seems to be relevant to me are slow modulations of the amplitudes, which could indicate some kind of crossfrequency coupling between different modes/processes.
Line 308: What do you mean by “response functions of the various components”?
Lines 336343: This is not a result of the present work and therefore completely irrelevant for the conclusions.
In general, the conclusions section leaves the same impression as the abstract, being full of technical details making the reader lost with trying to identify the key findings of the paper. In any case, this part reads very repetitive.
Can the differences between the Fourier analysis of Laskar et al. and the iSSA results for the same time series reported in this work be simply explained by the fact that one method explicitly assumes harmonic oscillations while the other does not?
Line 371: “Eccentricity has the main influence on insolation.” I would strongly doubt that the physical effect (in W/m2) of eccentricity does actually exceed those of the other orbital components.
Lines 384386: I think it is statistically meaningless to discuss apparent cyclic components with long periods that only fit two or three times into the total length of the studied time series.
Fernando Lopes et al.
Fernando Lopes et al.
Viewed
HTML  XML  Total  BibTeX  EndNote  

644  206  14  864  9  7 
 HTML: 644
 PDF: 206
 XML: 14
 Total: 864
 BibTeX: 9
 EndNote: 7
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1