Thank you for these useful comments. The responses to your comments are below.

- Line 160, data assimilation: Franke et al., 2020 is a good suggestion. We will cite it.

- Figure 1: The figure is showing the calibrated data from the Temperature 12k database.  Each record is represented as a colored horizontal line, with records arranged from north to south. The right y-axis helps show the latitude of the records. I will add additional text to the caption to clarify the right y-axis. As for Southern Hemisphere changes, the Southern Hemisphere mid-latitudes arguably do show a later timing of max temperatures than the northern mid-latitudes (compare the location of the black vertical lines), but it is noisy. This comparison would be aided by a larger network of Southern Hemisphere records.

- Table 1: Most of the proxy records come from the Northern Hemisphere, where winter insolation increased through the late Holocene and summer insolation decreased.  Because of these insolation trends, it is not surprising that more of the annual and summer records show cooling than the winter records.  It’s possible, however, that other changes would be seen in a more spatially-complete proxy network.

- Line 271: We will emphasize the decrease in boreal summer insolation. For reference, boreal summer insolation is typically ~450 W/m2 and was ~25 W/m2 greater in the early Holocene, a difference around 5%.

- Lines 326-327: We agree. We will reinforce this point in the sentence.

- Lines 340-342 and lines 535-545: That’s a good point: differences in proxy types may help explain some of the spatial variability in the Temperature 12k proxy database.  Regardless, Fig. 4 in Kaufman et al. 2020 also shows a fair amount of similarity between composites using different proxy types.

- Line 374: We will emphasize this.

- Lines 411-412 and section 4.3: The issue of spatial and seasonal biases in the proxy network is an important one. As you mention, Southern Hemisphere and winter records are under-represented in the Temp12k database, and it would be great to have more.  With that said, some reconstruction methods account for these biases: in the data assimilation, spatial covariance patterns are used to help infer Southern Hemisphere temperatures; in Kaufman et al. 2020 (“Holocene global mean surface temperature, a multi-method reconstruction approach”), proxies are composited into latitudinal bands and averaged together using the relative area of each band, which gives the Southern Hemisphere signal the same weight as the Northern Hemisphere signal.  In the data assimilation, covariance patterns are also used to account for proxy seasonality.  More proxies would certainly be helpful, however.  We will add additional discussion of biases to the paper.

- Figure 10: Yes, the variety of proxy temperature signals at 0-1 ka emphasizes the large amount of spatial variability in the proxy database.

- Line 463: Thank you. We will check to see whether additional new papers should be cited. Matero et al., 2017 is already cited.

- Citation Osman et al.: We will fix this citation.

Thank you for these comments. We’ve replied to your comments in order below.

1) Proxy records in the Temp12k database have a large range of temporal resolutions. Some ice core, speleothem, and two wood records have high resolution, although these archive types are in the minority. Additionally, a considerable amount of marine sediment, lake sediment, and peat records have mean resolutions lower than 200 years. While bioturbation may smooth some of this signal, our approach attempts to retain higher resolution information rather than smoothing the proxy signals up front. This results in a climate reconstruction which, while nominally decadal, simply represents the information contained in the proxy database. We will clarify this in the paper and stress that the decadal nature of the reconstruction does not imply that it contains robust decadal information. Instead, we want to use the data present in the database. This is similar to the approach used in Marcott et al, 2013 (“A Reconstruction of Regional and Global Temperature for the Past 11,300 Years”), which presents a temperature composite with a resolution of 20 years despite a median resolution of 120 years in their proxy data.

2) In our approach, the data assimilation is computed on a decadal timescale, so we need to have assimilated proxy values on the same timescale. To regrid proxy data to a decadal resolution, nearest neighbor interpolation was used to regrid data to annual resolution (e.g., the value for a given year is equal to the closest value), then binned to decadal. The intermediate step of interpolating to annual is meant to better account for sub-decadal data and decades mid-way between data points. A similar approach is used in some of the composites in Kaufman et al., 2020 (“Holocene global mean surface temperature, a multi-method reconstruction approach”). The approach is simply meant to represent multi-decadal data on a decadal timescale. With this said, we may take a different approach that also addresses age uncertainty concerns; see #3 below for more explanation.

To investigate how our interpolation approach affects the reconstruction, two reconstructions were compared. Both use 200 year bins. In one case, proxy data was interpolated between bins using nearest neighbor interpolation (this reconstruction is shown in Fig. B2h of the paper). In the other reconstruction, no interpolation was done. Both approaches lead to similar global-mean temperatures. The uncertainty range of the reconstruction using no interpolation was larger, but only by a small amount. The main difference between the reconstructions occurs near 11 ka, where many of the proxy records are cut off. The non-interpolated reconstruction has a rapid change in temperature here, which is likely not a climate signal, while the interpolated version is smoother.

3) Currently, age model uncertainty is not accounted for in the reconstruction, as mentioned in lines 542-543 & 586 of the paper. To account for age model uncertainties, we plan to rerun the reconstructions using a new methodology, which will work as follows: for each proxy record, we will generate a set of possible proxy realizations using the geoChronR R package (https://nickmckay.github.io/GeoChronR/) that sample the proxy’s age and magnitude uncertainties. Each proxy’s ensemble will be used to quantify the median proxy estimate and the joint age/magnitude uncertainty at a decadal timescale. (This approach will be used instead of the nearest neighbor interpolation approach described in the last comment.) The median proxy estimate will be assimilated and the joint age/magnitude uncertainty will be used to define a time-varying uncertainty term at each decade for each proxy (i.e., the R term in Eq. 2 of the paper). By combining age uncertainty and magnitude uncertainty into the R term, this method will account for age uncertainty without having to run a large collection of reconstructions. To our knowledge, this implementation hasn’t been used yet for paleoclimate DA, but should account for proxy uncertainties with a reduced computational expense. If there is a good reason why this is not a good solution, please let us know.

4) In the multi-timescale DA approach, the model data is used to compute covariances between different temporal resolutions. To use an example, let’s say that a proxy data point represents a 50 year mean. To relate 50 year means to the decade being reconstructed, the model-based proxy estimate is computed as a set of 50 year means, with each 50 year mean centered (to the extent possible) on the decades in the prior using a running average. In theory, this lets us quantify how 50 year means relate to decadal climate, and this information is used in the data assimilation. Additionally, since a data point representing a 50 year mean spans 5 decades, the process is repeated for each of the decades, with the only difference being our use of slightly different model values due to the moving prior window. The use of a running mean for multi-decadal averages is not perfect, but difficult to avoid without having a very large number of modeled years with which to construct the prior. In the paper, we will clarify the description of the multi-timescale data assimilation approach.

This multi-timescale approach is theoretically useful, as it is designed to allow each proxy to inform the reconstruction on its own timescale. While the final result is not drastically different from a reconstruction made using a single-timescale binned approach, this paper is partly an exploration of new methodological choices. We feel like that is a good justification for using the method and alternate experimental designs are shown in Appendix B. Regarding model covariances on different timescales, it’s possible that the models do not properly quantify climate covariances on different timescales; however, that would be a shortcoming of the data rather than the method. In the released code, if users do not want to use the multi-timescale approach, it can be disabled by setting two variables (time_resolution and maximum_resolution) equal to each other in the config file.

5) You’re correct that 0.09C may be too small to be considered “warmth,” as suggested by the uncertainty bands in Fig. 12, which include 0 from ~8 ka to the present. Using a t-test of two related samples of 1002 values each (the size of the prior), the mid-Holocene period (5.5 - 6.5 ka) is significantly different from the recent period (0-1 ka). However, it’s unclear if this is the right significance test to use. If ensemble members are analyzed individually, 88% of them have a positive anomaly at the mid-Holocene and 12% have a negative anomaly. We welcome better suggestions for a significance test. If a clear significance test cannot be established, the paper text can be modified to remove the assertion of mid-Holocene warmth, instead stating that our reconstruction does not support a cooler mid-Holocene.

6) For this first paper, we use the calibrated data from the Temperature 12k database, which is readily available for a wide range of proxy types. PSMs provide many potential advantages for data assimilation and represent an exciting avenue for exploration, but are not available for all of the proxies we use (e.g. pollon). Additionally, some PSMs are best used with isotope-enabled model simulations, which limits the choices of prior for the Holocene. The use of PSMs will be a focus of future study.

Line 65: We will be more specific about the inclusion of land records. However, it is difficult to describe the full difference in the databases, since Osman also includes many marine sediment records which are not currently in the Temp12k database.

Line 80: We will clarify the time-resolution and age control proxy criteria.

Line 98: We will include a figure of proxy temporal resolution.

Line 99: It would be useful to have a better sense of the temporal coverage of data points. However, that data is not widely available in the Temp12k database. As such, we’ve decided to treat data as continuous and mentioned our justification in the paper. While unrealistic, this assumption has the benefit of avoiding possible non-climatic signals generated by the alternation between proxy data and data gaps. To explore the effect of interpolating proxy data, we did the comparison described in point #2 above; however, we plan to replace the interpolation approach with the ensemble approach described in point #3 above.

Line 126: The minor checkerboard pattern is present in the results of the HadCM3 simulation. It seems to result from processes involving the ocean streamfunction and the model orography, although the details of this are not completely clear. In time-slice simulations, rerunning a simulation with smoothed orography causes the checkerboard pattern to disappear but the general climate remains the same. Because of this, we decided to simply apply a smoothing filter to the model results, which is mentioned in line 126 of the paper. We do not expect this to have much effect on the reconstruction.

Line 137: Yes, decadal resolution is higher than the mean resolution of most of the records. This is done by design, because we want to utilize information from our high-resolution records. Making use of data at a variety of temporal resolutions is one of the primary goals of the multi-timescale data assimilation methodology. As mentioned in point #1 above, the widely-cited Marcott 2013 proxy composite also uses a temporal resolution at the high end of the utilized proxy data. Regardless, we will add a note that our reconstruction is “nominally” decadal, but actually represents the information content in the proxies themselves.

Line 158: We agree this wording was unfortunate, since of course there are many PSMs that have been published by both you and by Sylvia and others. We will rephrase this statement. This first study uses the calibrated temperature records and we plan to move forward with exploration of proxy system models in the future. One complication to the use of proxy system models is that many use isotopes as inputs, and there is a limited availability of isotope-enabled simulations relevant to the full Holocene. Osman et al. 2021 used newly-run isotope-enabled simulations for their prior, but running new simulations is a considerable undertaking. Instead, we plan to explore the use of existing simulations, such as iTRACE. Additionally, by using calibrated temperatures for this first reconstruction, we can use it as a baseline to explore the added value of proxy system models in the future.

Line 515: This calculation uses the time windows of 5.5 - 6.5 ka vs. 0 - 1 ka. This is specified on line 515.

Line 518: This is discussed in comment #5 above.

Line 537: We will cite Osman et al., 2021 in this discussion.

Line 576: Unfortunately, metadata about sample thickness is not present in the Temp12k database for the vast majority of records. It would be good to have this information standardized and included in future databases.

Regarding localization: Covariance localization presents some benefits and some drawbacks, which are discussed in Appendix B.2. Among the drawbacks: 1) it is not clear which long-distance relationships are valid, so covariance localization is largely arbitrary and 2) covariance localization diminishes the overall impact of proxies on the reconstruction, allowing the prior to have greater influence on the reconstruction. However, the approach also has benefits, so we show some experiments using covariance localization in Fig. B2 and localization is an option in our released code. We will continue to explore its use in future reconstructions, and will clarify the discussion in Appendix B.2 if needed.

Line 835: Interesting! I didn’t realize that the covariance localization was possible with the simultaneous data assimilation approach. We will fix this oversight in the text.

Thank you for the review.  We respond to your comments below.

1. Our reconstruction differs from the one presented in Osman et al., 2021 in a variety of ways: the proxy selection, the use of PSMs, and the chosen model prior. Additionally, while we used some of the same proxy records as in Osman et al., 2021, they used additional marine records. The marine sediment records from Jess Tierney that we considered for assimilation were 103 d18O, 63 alkenone, 51 Mg/Ca, and 13 GDGT records, but they have added additional records to their database. We will look into acquiring more records for future work.

With this said, we ran two new experiments to explore the effects of land and ocean proxies on the data assimilation: one experiment only uses land proxy records and the other only uses ocean proxy records. These two experiments share some similarities, but differ in their spatial patterns especially in the early Holocene, where the ocean-only reconstruction is warmer. For the global mean, the default all-proxy experiment has a mean anomaly of 0.09°C for the 6-0.5 ka millenniums (as discussed in the paper) while the anomalies are 0.08°C and 0.04°C for land-only and ocean-only experiments, respectively. The low value of the ocean-only reconstruction suggests that some of the mid-Holocene warmth signal is coming from the land proxies, which are absent from Osman et al., 2021.

2. More Southern Hemisphere records would always be welcome, but data assimilation accounts for spatial biases to some extent. Essentially, data in unknown regions is estimated based on known data and model covariances. However, it’s likely that the reconstruction would be improved if it had more Southern Hemisphere records to rely on. The localization radius experiments in Fig. B2a-d are somewhat relevant to this question; in the default experiment, the large amount of proxies in the Northern Hemisphere may have an outsized effect on the reconstruction (after being translated through the prior’s covariance pattern), but in experiments with a localization radius, proxies cannot influence the reconstruction in distant regions. The global mean temperature reconstruction is not dramatically different between these experiments.

3. There are a few proxies in the network that do have published PSMs, notably ice cores, speleothems, and some of the tree data. However, we decided to use calibrated records for this first paper on the project, and move forward to a sensitivity test using the nonlinear/water-isotope-enabled PSMs for a second study. Furthermore, most of the proxies we could use PSMs for depend on water isotope enabled model output. The availability of water isotope enabled output is also very limited for the Holocene (e.g., iTRACE results have only been released prior to 11 ka and iCESM covers 850-2005), but nevertheless we could use the early Holocene and last2k isotope enabled simulations to conduct sensitivity experiments. The use of such PSMs will be a focus of future work in this project.  We will add additional discussion about this to the paper.

4. The proxy records initially have a variety of temporal resolutions, from sub-annual to multi-centennial. Since we are doing the data assimilation at decadal resolution, we need to have the proxy data at a decadal resolution. As mentioned in the paper, we make the assumption that all proxy data is continuous.  To represent longer-time-averaged data in a decadal context, we use a nearest-neighbor approach.  Data is first interpolated to an annual resolution using nearest-neighbor interpolation, then binned to decadal.  This two-step process is similar to doing a decadal nearest neighbor interpolation, but is meant to better account for sub-decadal data and time intervals mid-way between data points.

If proxy data was instead assumed to represent non-continuous periods (i.e., representing shorter time means with gaps of no data in between), the variability in the reconstruction may be higher, but some of this may be non-climatic variability resulting from records alternating between having data and lacking data. To test this, we ran two new data assimilation experiments.  For simplicity, both use 200 year bins. In the first case, data was interpolated. In the second case, no interpolation was done. In both cases, the global-mean temperature reconstruction looked similar. The main difference occurred near 11 ka, where many proxy records end; in the reconstruction without interpolation, there is greater variability here, but we do not think that such variability represents a real climate signal.

5. The two models used in the prior are on different grids, so interpolation is needed to put them on the same grid.  The chosen 64 x 96 lat/lon grid is somewhat arbitrary, but is in-between the resolutions of the two models (73 x 96 for HadCM3 and 48 x 96 for TraCE).

6. There are multiple reasons why one would use a multi-scale approach and as you say, one reason would be to capture longer term ocean memory or dynamics. For our reconstruction using multiscale assimilation is particularly important because we have proxy data on very different timescales, from multi-centennial sediment records to high resolution ice core records. Given timescales that differ by orders of magnitude, we wanted to treat each proxy on its own timescale and not assume that all proxies provide meaningful information on only a single timescale. It is a straightforward exercise to show that averaging or interpolating high resolution data down to low resolution removes high frequency variability and spectral information; this is precisely what happens when we assimilate proxy data that has been bin-averaged: high-frequency information is lost in the reconstruction. Additionally, there will be differences in the covariance structures in the prior at decadal vs. multi-centennial scales.  We will revise the main text to discuss this issue more fully.

7. We will add additional explanation of how the data assimilation methodology handles seasonal proxy records.

8. As mentioned in comment 1 above, Osman includes marine records which were not included in our reconstruction. Also, because the two methods differ in their prior and use of PSMs, a straightforward comparison is not possible.  The comparison between land-only and ocean-only proxies described in comment 1 above partly addresses this topic.

9. To look into this, we did a test similar to the one on Extended Data Table 1 of Tierney et al., 2020 (“Glacial cooling and climate sensitivity revisited”). A subset of proxies was withheld from the data assimilation, then correlation, CE, and RMSE values were calculated for un-assimilated proxies at each age. As in Tierney et al., proxy scaling factors of 1 (the default), 1/2, 1/5, 1/10, 1/20, 1/100, and 1/500 were tested. Median verification metrics were best for values of 1, 1/2, and 1/5. Since none of the other scaling factors were unequivocally better than the default, we will continue to use the default values.  However, a reconstruction using a 1/5 scaling factor is still shown in Fig. B2 of the paper.  As for ECR and ECPSS, we are unfamiliar with those methods. If you provide further explanation for how those methods would help, we can investigate their use.

10. Reconstructing additional climate variables will be a focus of future work.  Before providing reconstructions of those variables, we would like to assimilate more proxy records.  In particular, we are currently working to assimilate hydroclimate proxies.