General comments:

The paper presents an interesting analysis of selected records and is highly relevant to the study of centennial to millennial scale natural variability, providing some valuable insights into the space-time dynamics of the North American climate over the Holocene. I think that most of the prominent features discussed are convincing, but I think that limitations of the data and analysis can explain why little coherent variability is observed otherwise.

First, I think that a more formal definition of cen-mil variations could be given. When amplitudes of variations are given, it should be clear that those are (if I understood correctly) the standard deviation of the 50-year interpolated series over the last 7000 years after removing a loess trend with span of 6/11. If that is indeed the definition usually assumed, then I would question it since lots (or virtually all) of the records will have a resolution lower than 50 years, and the resulting standard variation will be underestimated. In addition, time-uncertainty will further smooth variability when taking such averages after interpolation. Therefore, the amplitude of the resulting cen-mil variations can only be seen as a lower bound for true variations (assuming that non-climate variations are small and mostly averaged out, which might not be exactly true, but I think this is still very challenging to quantify).

Quantifying variability at a given timescale could also be achieved by integrating the power spectrum over a given timescale/frequency band to obtain the variance. That said, removing the LOESS is an effective high-pass filter (although the frequency cutoff is not precisely defined) and will give similar results, with the advantage of being in temperature units. Having variability in temperature unis could also be achieved using Haar fluctuations (Lovejoy & Schertzer, 2013), which incidentally works quite well with irregular data (Hébert et al., 2022, NPG).

Second, the null expectation model could be improved. On line 156-157, it is stated that LOESS filtering white noise gives time series with temporal autocorrelation characteristics similar to paleoclimate timeseries, well, not really. The LOESS is quite an abrupt low-pass filter, so if you apply it to white noise, the result is mostly white noise at a lower resolution with a smoother interpolation between, the temporal autocorrelation isn’t changed above a certain timescale. Why not generate data that already has realistic autocorrelation (e.g. fractional noise with an exponent estimated from paleo data, which often hovers around beta=1/H=0) and then low-pass filter it (as in Hébert et al., 2021, NPG for example)? I would consult Reschke et al. (2019): https://doi.org/10.1016/j.cageo.2018.11.009, which wholly focus on the method aspect of creating an appropriate null expectation model and provides tools for testing the significance of correlations when dealing with irregular paleoclimate data.

I think the null expectation model used by the author is pedagogically useful and helps getting a conceptual understanding of the characteristics of paleoclimate data. It is not however a state-of-the-art representation of paleoclimate data, I would say.

Detailed comments:

Section numbers 1.1 and 1.2 should read as 2.1 and 2.2 I believe.

Figure 2: The shading should be described, are we looking at the spread between the single curves used in the average?

Line 254: Unclear to me where the range of random expectations comes from, Fig. 3a, which is cited, is not showing that. Are we talking about the histogram on Fig. 1a? Or are we averaging a different white noise on both series and calculating the correlation 1000 times and using this as the range as in Fig 1b? I would consider using the “corit” package from Reschke et al. (2019) to calculate the correlations and their significance.

Line 264: Again would be important to specify how the range is obtained.

figure 3B: The blue curve appears to go outside of the shading around 1000 year BP or so. Unclear how the shadding is calculated.

Line 286: Meaning of sentence is unclear to me.

Figure 4: I would write the correlation coefficient with the same size font as the information on the strength of the correlation is already in the value itself. Maybe adding a colour code for whether they are significant or not would be useful. I think the figure could generally be improved. Are we learning a lot from the scatter plots and the distribution of detrended residuals?

Figure 5: The use of shading instead of a curve to reproduce the MAT spectrum is confusing I'd say, why not keep the black curve reproduced? There is likely very little signal left on these curves since we are looking at a temporal average of  curves with time-uncertainty (see e.g. Hébert et al., 2022, Nat Geosc, Extended Data Fig. 1 for a surrogate experiment with time-uncertainty), which smooths out drastically the high-frequency signal, and the LOESS detrending removes the low-frequency signal. The steep scaling on the high-frequency side is very likely the result of power-loss due to combining time-uncertain series. I would only expect a narrow band between 2000-4000 (maybe 1000-4000) to be somewhat reliable. Also, power spectral density is not unitless, I would advocate showing the actual units rather than the log (but keeping log spacing); as a result, it is unclear to me how the series with different units are compared, was the spectrum computed on something like the z-score of the residual timeseries?

One last point on this figure, the uncertainty given for the scaling exponents are quite meaningless, and I would advocate giving approximate values (e.g. beta≈0.7). There is nothing wrong with taking a linear fit of logs to get an approximation of the scaling exponent, but many assumptions that do not hold are made when doing so, e.g. about the distribution of the residuals being independent and normal. If we really wanted to make meaningful estimates of uncertainty here, I would recommend designing a surrogate model for the series and using a Monte-Carlo approach (e.g. Hébert et al., 2021, NPG), ideally also accounting for time-uncertainty and other proxy processes (Kunz et al., 2022, Clim Past; Dolman et al., 2021, Clim Past). However, I would also be satisfied by approximate values without any uncertainty given (e.g. beta≈0.7).

Figure 6: B) The map could use a colour scale. I would also use only two colours on a divergent scale (blue to red for example) rather than orange-green-white-blue-purple. The colour scale for the colour of the triangles would also be useful. A) The shading should be described, are we looking at the spread between the single curves used in the average? In the caption, I would quickly mention those are the ENA sites.

Line 358: Not sure I am convinced by the significance of the "Evidence of the 8200 YBP cool event" on Fig. 6A.

Line 359: Maybe "warming phase" would be more accurate than "warm phase”, although it is not warming much, but rather stays constant until it rejoins de long-term cooling trend. I would generally avoid describing minor wiggles as I am not convinced they are robust.

Line 364: How robust is this correlation? It seems to me like there is very little correlation outside the first 1000 years (when we have a very steep change) and the 5500-4800 period, which was used to discriminate the records. Could the apparently event in the 5500-4800 period be produced by construction? Say we would take the same set of records and average them based on their trends in another arbitrary time period, would’t we also get a similar result but for a different time period, i.e. the records with (randomly) positive trends averaged together will create an average series with a stronger trend over the selected time period and average out elsewhere? I did a quick experiment and simulated 20 fractional noise series with beta=1/H=0 (following out method in Hébert et al., 2021, NPG) covering 11ka, low-pass filtered and sub-sampled at 250 years resolution, then calculated trends over a 1000-year long period in the middle the series, and then stacked the series with a positive trend together, and the ones with negative trends together. This yielded an average correlation of -0.14+/-0.27 between the two stacks, and -0.15+/-0.17 between the two stacks after a 6k loess trends was removed in each. Thus, I expect that the method may induce a negative bias on the correlation and would suggest to run a similar surrogate experiment to assess the significance of the correlation, and also to try out stacking the records with respect to the trends estimated over different reference period and see whether we can easily get such events.

Line 367: How do we know that this standard deviation difference is due to different variability? Then why are the warming/cooling events around 5500 BP of similar magnitude? Isn’t it possible that they are the result of different different resolutions (and also time-uncertainty) in the input timeseries? It would be easy to check whether the mean resolution of the northern records is lower than that of the southern ones.

Line 369: Why is this slope relevant?

Line 384: “appear evident” rather than “appear evidence”? Revise sentence.

Line 385: Unclear what are those small temperature fluctuations that are discussed, is that with respect to the standard deviation of individual detrended series after interpolation at a 50-year resolution? Seems to contradict line 239 where 0.25-0.5 variations are discussed. Cen-mil variations could be formally defined.

Line 391: Not sure that the null model can be described as autoregressive, which is defined (on Wikipedia) by: “ The autoregressive model specifies that the output variable depends linearly on its own previous values and on a stochastic term (an imperfectly predictable term); thus the model is in the form of a stochastic difference equation (or recurrence relation which should not be confused with differential equation).” As the LOESS filter is not causal, this means that it also depends on future values.

Line 409: I would see this drop in high-frequency variability as a consequence of averaging time-uncertain records together, having bioturbation or sediment mixing, and interpolation of lower resolution sediments.

Figure 8: MBM is written mbm on the axis of the top panel.

Line 501: Not sure if this statement is substantiated.

Line 504: I’d say “relatively” well-dated.

References:

Dolman, A.M., Kunz, T., Groeneveld, J., Laepple, T., 2020. Estimating the timescale-dependent uncertainty of paleoclimate records – a spectral approach. Part II: Application and interpretation. Climate of the Past Discussions 1–22. https://doi.org/10.5194/cp-2019-153

Hébert, R., Herzschuh, U., Laepple, T., 2022. Millennial-scale climate variability over land overprinted by ocean temperature fluctuations. Nat. Geosci. 15, 899–905. https://doi.org/10.1038/s41561-022-01056-4

Hébert, R., Rehfeld, K., Laepple, T., 2021. Comparing estimation techniques for timescale-dependent scaling of climate variability in paleoclimate time series. Nonlinear Processes in Geophysics Discussions 1–26. https://doi.org/10.5194/npg-2021-7

Kunz, T., Dolman, A.M., Laepple, T., 2020. A spectral approach to estimating the timescale-dependent uncertainty of paleoclimate records – Part 1: Theoretical concept. Climate of the Past 16, 1469–1492. https://doi.org/10.5194/cp-16-1469-2020

Lovejoy, S., Schertzer, D., 2012. Haar wavelets, fluctuations and structure functions: convenient choices for geophysics. Nonlinear Processes in Geophysics 19, 513–527. https://doi.org/10.5194/npg-19-513-2012

Reschke, M., Kunz, T., Laepple, T., 2019. Comparing methods for analysing time scale dependent correlations in irregularly sampled time series data. Computers & Geosciences 123, 65–72. https://doi.org/10.1016/j.cageo.2018.11.009