General comments

We would like to thank the reviewer for the comments and provide our responses and proposed changes to the manuscript in italics below:

Peatlands are important carbon reservoirs at the Earth’s surface; their shrinking or thickening, partially affected by human activities, can contribute to carbon release or sequestration and thus exert an impact on climate to some extent. This manuscript by Moore et al. focuses on the evolution of upland peat thickness in a slowly eroding area in Belgium, using paired 10Be-26Al in underlying saprolite and nearby stream sediments. The 10Be-26Al concentrations suggested lower values compared to steady-state conditions at current peat thickness, hinting at a thicker cover in the past. The authors attributed the shrinkage of the peatland to enhanced human activities. Note that I am not an expert in peatlands or burial dating, and thus I may not be able to comment on part of the paper. In general, the utilization of cosmogenic nuclides in tracking peat thickness based on its shielding effect on underlying saprolite seems to be a novel application. Nevertheless, such an application relies on some assumptions, as stated in the manuscript, and the degree of uncertainty they may introduce is difficult to quantify, which might affect the robustness of the data interpretation. I hope the authors find the comments below useful during their revision process.

Specific comments

1. Does this application rely on an assumption of time-invariant saprolite denudation? However, if peat thickness can vary—for example, from being completely removed to a thick cover—the denudation rate may also vary to a large extent, ranging from fast erosion by runoff to slow erosion due to the shielding effect of the peat cover?

2. The study area is characterized by a MAT of 6.7 °C. As such, I would assume that the temperature in the LGM could be very close to the freezing point (mentioned in the manuscript). In this case, even if there was no glacial coverage (please provide literature) in the past (e.g., 1 Myr), periglacial processes could play an important role in surface denudation, particularly during past glacial periods. How does this process potentially affect your modeled denudation rate?

This application does not rely on saprolite denudation rates being time-invariant, only that changes in denudation rate occur with a periodicity that is short in comparison with the integration timescale of the signal. The integration timescales range from 0.7 to 1.8 Ma (Table 3), whereas changes in denudation rate are most likely to occur over glacial-interglacial (100 ka) (in the case of periglacial processes) or shorter timescales. Thus, changes in nuclide concentrations from climatically-driven changes in denudation rate should be highly damped and the rate we infer is essentially a time-averaged rate.

We will add a sentence describing this to section 4.3 where we explain the integration timescale.

3. Eq. (2). The calculation of the long-term denudation rate using water chemistry is indeed inappropriate. Beyond the timescale bias and effect of mineralogical sorting mentioned by the authors, the solute sources can be derived from peat, quartzite, and argillite, while the solid source only includes one source (either quartzite or argillite). This discrepancy might make the application of Eq. (2) for denudation rate calculation problematic.

This is a good point, although it seems that further weathering of primary minerals in the peat is probably low. Dry peat generally has a density of around 0.1 g cm^-3 and is about 10% mineral ash by weight. Thus, even meters thick peat corresponds to only a few g per cm^-2 of mineral material that could weather and release Si.

Even if Si release from the peat were significant, this would mean that we have over-estimated the denudation rate of the saprolite by attributing all the dissolved Si flux to saprolite weathering. A lower denudation rate of the saprolite would be in even better agreement with the denudation rates inferred from the cosmogenic nuclides.

In the revised manuscript, we will try to frame the analysis of denudation rate from water chemistry and solid-phase mass balance more explicitly as an order-of-magnitude check on the cosmogenic nuclide results rather than a central line of evidence.

4. Line 188. Is there a seasonal change in water content, and how does it affect your estimate of atomic mass?

There is some seasonal change in water content at the site but it is quite modest (Henrion et al., 2025). The peat likely remains saturated below ca. 30 cm depth except under extreme drought conditions.

Even in totally desiccated peat, hydrogen is still abundant because it is a major component of organic matter. Thus, the mean atomic mass of dry peat is likely closer to water than rock.

We will add a sentence in section 3.5.1 describing this.

5. Line 292: The large difference in immobile element ratios between saprolite and bedrock is indeed worrisome. But I agree it is difficult to find a convincing explanation.

6. Conclusions: Lines 577-578. Since the authors argue that the removal of peat is mainly caused by human land use in the last few hundred years, why not use a shorter-half-life nuclide such as 137Cs (Jelinski et al., 2019, JGR), which may be more appropriate for this purpose? Note that I do not mean to suggest measuring an additional nuclide, but rather to provide a short paragraph on future directions.

It would certainly be interesting to consider the dynamics of shorter-lived fallout nuclides in this system. However, this is beyond the scope of this work, as it requires not only additional measurements, but different conceptual and mathematical models. We prefer to keep the message of the paper focused on applications of in situ cosmogenic nuclides in peatlands, which, to the best of our knowledge, have not been previously explored, but we do mention in situ ¹⁴C as another interesting direction to explore. Moreover, bomb-pulse ¹³⁷Cs is limited to studying processes since the mid-20th century. We believe the land-use history at our site extends further back in time and would thus be difficult to fully resolve using bomb-pulse nuclides.

Technical corrections

7. Fig. 1: Perhaps add a colorbar of elevation for the inset map of Belgium. 8. Table 1: Depth to top or bottom. What does top or bottom mean? Top of the peat or top of the saprolite?

We will add a color bar and clarify top and bottom by adding a note in the table caption.

We would like to thank Richard Ott for the helpful review. If the manuscript is accepted for revision, we would be pleased to add clarifications and expanded justifications of the inversion approach and to implement the other suggestions.

Moore et al. present an interesting study trying to infer past peat thickness, denudation, and timing of peat accumulation in the Belgian Ardennes. The study is well written and in most parts also very thorough. However, parts of the methods were brief and thus I could not follow exactly what the authors did to run their MCMC inversion. Based on my understanding of what the authors did, I have concerns that the recovered parameters  and the conclusions might be biased.

Inversion of CRN concentrations:

My major concern is the application of the MCMC inversion for this study. The methods section on the MCMC application was brief, and thus, I might have misunderstood something. It was not clear what parameters were inverted for, what the prior ranges were, and whether samples were inverted jointly or separately. The way I understood the manuscript, the authors use the MCMC to constrain the parameters of denudation rate, initial overburden thickness, overburden thickness, and time of change, in equation 5 INDIVIDUALLY for every hillslope sample. Thus, there are 4 unknowns for every sample (D, z1, z2, t), but only two equations with data (10Be equation and 26Al equation). This is an undetermined problem and cannot have a unique solution without additional information.

This can be somewhat seen in the lower panels of figure S1, where one would expect all parameter combinations that can explain the data to fall onto a curve. They do so for panel D , but the MCMC apparently did not explore the entire space for denudation rates. From personal experience I know that these MCMC random walks sometimes tend to cluster one side of the prior parameter space—for reasons I don’t understand—, leading to the effect that they don’t map the full line of parameter combinations. However, in theory, with 3 parameters one can calculate curves of valid parameter combinations analytically from the two nuclide build-up equations.

The only reason to apply a MCMC to this problem would be to invert all samples together, hoping that differences in the sample parameter combinations provide additional information. This is the approach that Hippe et al., 2021 used for 3 parameters and 2 nuclide equations. However, they unfortunately did not provide synthetic tests to prove that the inversion can still reliably recover synthetic parameters (I think it’s great that the authors of this study included such a test!). My experiences from personally working with similar inversions for cosmogenic nuclide erosion histories is that the inversion only partially recovers initial conditions, for combined samples with 2 nuclides and 3 parameters.

A simple fix could be to hold values constant—assume a denudation rate and/or zero initial peat thickness—and then calculate the other two parameters for every sample, resulting in a range of plausible parameters.

This being said, I apologize, if I misunderstood part of the manuscript and the samples were indeed jointly inverted. In that case, some method clarification and a better documented synthetic text with several samples would be nice.

Introduction

We appreciate the critical evaluation of the MCMC inversion and agree that this is an under-determined problem since we invert for 4 parameters from 2 nuclides. It is true, a unique solution does not exist. Yet, one can still determine the joint probability distribution (i.e., the probability that certain parameter combinations occur together) using MCMC.

Our samples are indeed inverted individually rather than jointly. This is because our model has no parameters that are constant between samples (i.e., we allow there to be a unique denudation rate, time of change, and different overburden thicknesses before and after the change for each sample). This differs from the main approach in Hippe et al. (2021) where the timing of the change in erosion rate is assumed to be constant across all samples. However, an individual inversion where the erosion rate is allowed to vary between samples is also presented in their supplement (their Fig. S4).

Markov Chain Monte Carlo has been applied to find joint distributions for cosmogenic nuclide models with more parameters than samples in a number of studies from the last 10 years or so, but, as rightfully pointed out, is not without limitations and is sensitive to the choice of priors (we will move our priors from the Zenodo repository to the main text to make them more explicit). It is an interesting question – how well does the MCMC find the “true” parameter values for our models?

Least-squares gives similar results to MCMC for LS2-LS5

One way to evaluate the MCMC is to simplify the model to involve only 2 parameters and solve using classical least-squares. Because the MCMC results indicate that the initial overburden thickness and timing of a change in over-burden thickness are unconstrained, it makes sense to neglect these two terms. If we do this and solve by least-squares, then we get the following results for peat thickness (in g cm^-2):

and for denudation rates (in tons km^-2 yr^-1)

All results from the least-squares inversion fall within the MCMC interquartile range. Thus, the peat thicknesses and denudation rates determined from LS2-LS5 are not artifacts of the MCMC inversion.

Why use MCMC at all? The reason is that it is not possible to explain the concentrations measured in LS1, which falls below the secular equilibrium line on the two-isotope plot, using a 2 parameter steady-state exposure model. There must have been a change in shielding far enough in the past for decay to be important. Although we do not really need the MCMC for LS2-LS5, for consistency we use the same model and inversion method for all the hillslope positions. Likewise, the catchment-averaged sediment model contains more than 2 parameters, and so a unique solution is also not tenable there, causing recourse to MCMC. 

The case of LS1

Are the parameter estimates from the MCMC for LS1 reasonable, and why does Figure S.1c not show the expected trade-off curve between denudation rate and overburden thickness?

There are several lines of evidence indicating that the inversion results for LS1 are reasonable:

1. The extremely low denudation rate is consistent with the paleo-surface model from qualitative geomorphology.

2. The inferred long-term peat thickness from LS1 is in line with the other samples and gives a reasonable hillslope profile (Fig. 3). Because the samples were inverted individually, these are independent measurements.

3. We have independent knowledge that uplift occurred around the modeled time of peat thickening. This provides a plausible explanation for the thickening.

With respect to the denudation rate – overburden thickness trade-off curve (Fig. 1c), part of the explanation likely lies in the fact that the denudation rate cannot go too high. This is because as denudation rate increases the integration time decreases, and the signal becomes less able to capture the radioactive decay needed to generate ratios below secular equilibrium. Also, the values for the denudation rate are not controlled only by the trade-off with overburden thickness, but also interaction with the timing of the change and the initial overburden thickness.

Proposed changes to the manuscript

To clarify our approach we propose making the following changes to the manuscript:

1. We will move the priors from the data repository into a table in the manuscript. This will make more evident the values of the priors that we used and which parameters we are inverting for.

2. We will provide more discussion of our MCMC methods in the text (section 3.5.4) and state explicitly that the samples are inverted individually.

3. We will note that essentially the same results are obtainable for LS2-LS5 using classical least-squares methods. We will include the least-squares results in the supporting information and add the code used to calculate these to the Zenodo repository. We will also expand the caption to Figure S.1 to better describe the synthetic models.

Other comments:

Production rate scaling. The authors use Stone 2000 scaling. However, compared to the last few million years the geomagnetic field is comparatively strong at the moment. This may lead to an underestimation of long-term production rates. Maybe given the other uncertainties in the inversion this effect is small enough to be neglected. In that case, the authors could add a brief statement explaining their reasoning.

At the study site, which is located above 50° N, the latitude effect is muted. The difference between Stone (2000) scaling factors and the Balco et al. (2008) time-dependent factors averaged over the last 1 million years calculated using CRONUSCalc (https://cronus.cosmogenicnuclides.rocks/2.1/) seems to be <1% for the study site.

There is a larger offset between the production rate predicted by the Lifton et al. (2014) scaling model over the last 1 million years (ca. 6.6 atoms g^-1 yr^-1 for ¹⁰Be) and Stone (2000) (ca. 7.2 atoms g^-1 yr^-1).

However, the impact on the parameter estimates is quite small. This is because under peat thicknesses greater than 260 g cm^-2 (the lowest median modeled value) the high-energy nucleon flux is severely attenuated. This is especially true when considering a mean free path of 110 g cm^-2 through peat, which implies spallation is attenuated to <10% of the surface value under 260 g^-1 cm^-2. Muon-induced production therefore contributes disproportionately to total production.

Recalculating LS2-LS5 using the least-squares approach and a spallation production rate of 6.6 atoms g^-1 yr^-1 gives overburden thicknesses as follows (g cm^-2):

and for denudation rates (tons km^-2 yr^-1):

These results differ by only a few percent. We will add a sentence or two discussing the relative insensitivity to the spallation production rate to section 4.

The 10Be and 26Al measurement data need to be reported in more detail to ensure reproducibility. Currently, only the final concentrations are reported in table 1. However, the authors should report all relevant data that were used to convert AMS measurements to concentrations (AMS ratios, sample weights, carrier concentration, blank values, etc.).

All of this data (ratios, weights, spikes, blanks, Al concentrations, etc.) is reported in the supplementary data tables in the Zenodo repository. We will either add a statement in the table caption to make it clear that this is where the full data is, or else add another table to the manuscript, whichever is most consistent with ESurf’s policies.

I was confused as to why the authors only focus on Si for the weathering rates. Why weren’t the total dissolved solids considered? Why was the standard CDF equation modified to Si, given that supposedly all other major elements were measured by ICP-OES?

We chose to focus on Si because it is more likely derived from rock weathering than the other dissolved elements in the stream water. Stream solute concentrations at the site are low and cyclic salts dissolved in precipitation constitute a significant fraction of the total dissolved mass flux. For example, nearly all the Na⁺ in the stream water comes from precipitation, even though Na⁺ is the major element with the highest concentration in the water.

The standard way of dealing with cyclic salts is to assume that Cl^- is derived only from precipitation, and that the element:Cl^- ratio in precipitation is the same as seawater and correct accordingly. However, we lack the necessary Cl^- data to do this.

Si is the most abundant element in the bedrock, was present at a high-enough concentration in the stream water to be reliably determined by ICP-OES, and has an essentially negligible concentration in precipitation. Thus, to estimate the bedrock denudation rate without knowing the cyclic salt corrections, it is best to focus on Si. Of course, there is considerable uncertainty with this approach. Thus, we frame the solute-flux based denudation rates as a sort of “order of magnitude” check on the cosmogenic nuclide denudation rates. We will make this point more explicit in the text.

Would be nice to add a picture of a soil pit and sample location within the pit.

Attenuation length: The authors calculated a separate peat attenuation length. This makes a lot of sense. However, the peat depth according to the authors is only 0.2-2.1 m at the moment. Meaning that a lot of production will occur under the peat and thus have a different attenuation length. Does the forward model in the inversion use both attenuation lengths (for peat and below peat)?

Yes, the forward model includes this difference in attenuation length. We will add this point to the text in section 3.5.2.

L13: Seems like there’s a word missing after “explore”.

L122: What was the grain size of the stream sediment samples?

L159: Here or in the cosmogenic nuclide result table, report the carrier concentration.

L202: I suggest reporting the exponential parameters here or in the supplement. Otherwise, interested people would need to dig through the code.

L214: Missing i- for production pathway

L237: I presume the authors mean denudation rate, instead of erosion.

L440: Please, elaborate on your argument. The authors note the “nature of these lenses” make them unlikely to have had sufficient shielding effect. This is a vague statement that could be explained better.

Many thanks for reviewing the text in detail. We would be pleased to make all these corrections.

For questions about this review, feel free to get in contact.

Richard Ott