Bayesian inversion of a CRN depth profile to infer Quaternary erosion of the northwestern Campine Plateau (NE Belgium) 2

. The rate at which low-lying sandy areas in temperate regions, such as the Campine plateau (NE Belgium), have been eroding during the Quaternary is a matter of debate. Current knowledge on the average pace of landscape evolution in


Introduction
The Campine area is a sandy region which covers part of northeastern Belgium and the southern Netherlands (Fig. 1).It is part of the European sand belt and is drained by rivers that belong to the Scheldt Basin.The Campine area roughly coincides with the geological Campine Basin, the southeastern part of the North Sea basin.From a geodynamic point of view, the Campine Basin is located in an intermediate position in between the rapidly subsiding Roer Valley graben in the north, and the uplifting Brabant and Ardennes massifs in the south (Fig. 2).The Campine Basin has a long Cenozoic burial history.Post-Rupelian marine and estuarine deposition during the last 30 Myr almost exclusively consists of (glauconite-rich) sand, up to 300 m thick (Vandenberghe et al., 2004).From the Early to Middle Pleistocene onwards, Published by Copernicus Publications on behalf of the European Geosciences Union.terrestrial conditions became dominant with deposition of a thick series of fluvial sand and gravel from the Meuse and Rhine rivers (Figs. 2 and 3).In contrast to what the basinal setting of the Campine region would suggest, distinctive topographic features are preserved in the landscape.An illustrative example is the Campine Plateau, which shows a topographic relief of ca.50 m relative to the surrounding areas (Fig. 3).To date, quantitative data on the amount and rate of Quaternary erosion of the Campine landscape and the Scheldt Basin in general are missing.This stands in contrast to the availability of long-term erosion data from, for example, in situ-produced cosmogenic nuclides for the Meuse and Rhine basins (e.g., Schaller et al., 2001;Dehnert et al., 2011;Rixhon et al., 2011).Such data on catchment-wide erosion rates at multi-millennial timescales are crucial for determining background geological erosion rates to evaluate anthropogenic morphodynamics (Vanacker et al., 2007a), to provide calibration data for landscape evolution models (Bogaart and van Balen, 2000;Foster et al., 2015;Campforts et al., 2016) and to assess the overall stability of the landscape in the framework of long-term management of radioactive waste (Van Geet et al., 2012).
Cosmogenic radionuclides (CRNs) have proven useful for quantifying geomorphological processes over time spans covering the last 2 Myr (Schaller et al., 2001).Geomorphological surfaces can be dated by measuring the concentration of in situ-produced cosmogenic nuclides (e.g., 10 Be and 26 Al) that accumulated at the Earth's surface (Dunai, 2010;Hancock et al., 1999).As the observed cosmogenic nuclide concentration of a given outcrop is a function of its exposure age and denudation rate, stable (i.e., non-eroding) landforms provide optimal sampling locations for exposure dating (see, e.g., Rixhon et al., 2011).Most landforms are subject to erosion during exposure, resulting in a decrease in the cos-  eu/data-and-maps/data/wise-large-rivers-and-large-lakes) and location of the Scheldt Basin, (dashed line), the Nete catchment (dotted line) and the Flemish Valley (solid line).The general paleohydrography of the Meuse and Rhine between 0.5 and 1.0 Ma is shown in colored lines.Headward erosion as an explanation for the development of the Nete catchment is indicated with a yellow arrow.Beerten, 2005, andDeckers et al., 2014).mogenic nuclide concentrations with an increasing surface denudation rate (e.g., Dehnert et al., 2011).Braucher et al. (2009) showed that the exposure age (and post-depositional denudation rate) of eroding landforms can be constrained based on a deep (> 1.5 m) depth profile of a single cosmogenic nuclide that is sampled at regular intervals.

Fig. 4a
The accumulation of in situ-produced cosmogenic nuclides in eroding surfaces is a mathematical function with three parameters that are typically unknown a priori: the post-depositional denudation rate, E (m Myr −1 ), the exposure age, t (years), and the inherited concentration or inheritance, N inh , (atoms g −1 ).Unknown model parameters can be estimated by inverse modeling of CRN concentration vs. depth profiles.In this procedure, one iteratively proposes new parameter values until the model fits the observed data up to a given precision.This has been done for estimating E and t by, for example, Siame et al. (2004) and Braucher et al. (2009).However, model and measurement errors together with (measurement) data scarcity introduce considerable uncertainty in the optimized model parameters.The method by Braucher et al. (2009) accounts to some extent for analytical measurement errors as it generates several CRN concentration profiles consistent with the (analytical) measurement errors, computes for each model parameter set (e.g., a t − E pair) within a grid search the corresponding data misfits and retains the median misfit as the performance associated with a given parameter set.This allows one to derive a robust unique solution, but it does not quantify model parameter uncertainty and ignores model errors (that is, the model is assumed to be perfect).To assess model parameter uncertainty, Hidy et al. (2010) proposed performing plain Monte Carlo (MC) sampling from the pre-specified prior parameter distributions, ranking the resulting solutions according to fitting performance and retaining a certain percentage (typically 5 %) of the best performing solutions to compute parameter uncertainty estimates.
A more comprehensive quantification of parameter and prediction uncertainty is provided by the Bayesian framework.This approach uses Bayes' theorem to represent parameter uncertainty by a multivariate "posterior" probability distribution.The latter is given by the (normalized) product of a "prior" probability distribution, which represents available prior information, with a "likelihood" function, that encodes the deviations of the simulated (CRN) concentration data from the measured ones.Providing that the assumptions underlying the likelihood model are met, the posterior probability density function (PDF) contains all necessary information about the inferred parameters.Marrero et al. (2016) implemented Bayes' rule into the CRONUScalc program to derive the posterior parameter PDF.The approach taken by Marrero et al. (2016) is based on a MC variant where sampling is performed over a regular, 3-D lattice covering the prior ranges for E, t and N inh .It therefore requires a sufficient grid resolution to minimize the risk of missing E − t − N inh combinations with substantial posterior density.More importantly, the formulation by Marrero et al. (2016) considers solely the CRN measurement error(s) as source of uncertainty.This is theoretically valid only if the CRN model can fit the measurement data within the measurement error(s).In this work, we derive the posterior parameter PDF using state-of-the-art Markov chain Monte Carlo (MCMC) simulation (see, e.g., Robert and Casella, 2004), accounting not only for measurement errors but also (under certain assumptions) for model errors.Furthermore, we compare our Bayesian inversion approach with that of Marrero et al. (2016), illustrate the similarities and differences between the two approaches and propose a simple fix for making the uncertainty estimates from Marrero et al. (2016) more consistent in the presence of model errors.
The overall objective of this study is to infer within a Bayesian framework the potential post-depositional denudation of the northwestern Campine Plateau.This part of the Campine Plateau is drained by the Kleine Nete river, which belongs to the larger Scheldt Basin.It is an interesting test case because the northwestern edge of the Campine Plateau is covered by coarse gravelly unconsolidated sand from the Early-Middle Pleistocene Rhine and thus constitutes a fluvial terrace for which the depositional age nor the exposure age is well constrained (Beerten et al., 2017).

Geomorphological evolution of the Campine area
The post-marine hydrographical evolution of the Campine area started with the final retreat of the sea during the Neogene, as a result of systematic sea level lowering and overall uplift of the bordering areas around the southern North Sea (Miller et al., 2005;Cloething et al., 2007).During the Early Pleistocene, the Meuse followed an eastern course from Liège to the region north of Aachen where it merged with the Rhine (Fig. 2).Tectonic movements along Roer Valley graben faults, and uplift of the northern margins of the Ardennes-Eifel massif caused the Meuse to breach its northern interfluve and to follow a completely different course.At the same time, the Rhine shifted its course as well, flowing into the northern part of the Campine area where it merged with the Meuse (Fig. 2).Age control is limited, but this event probably took place around 1 Ma at the earliest, since the confluence area of both rivers was situated in the southeastern part of the Roer Valley graben prior to 1 Ma, and the area that covers the Campine Plateau today was drained by local "Belgian" rivers until that time (Westerhoff et al., 2008).It appears that both rivers shifted their course towards a more eastern position by 0.5 Ma at the latest, given the absence of Rhine deposits younger than 0.5 Ma in the depocenter of the Roer Valley graben (Schokker et al., 2005).The deposits that cover the Campine Plateau are often correlated with the upstream main terraces of the Meuse and high terraces of the Rhine (Paulissen, 1973).Westaway (2001) provides a time window for deposition of Rhine sediments west of the Ville Ridge (high terraces HT2 and HT3) between 0.5 and 1 Ma.The Rhine sediments on top of the Campine Plateau have been attributed to the Sterksel Formation, which was deposited between ca.0.6 and 1.1 Ma according to van Balen et al. (2000), and between ca.0.75 and 1 Ma (post-Jaramillo Early Pleistocene) according to Gullentops et al. (2001).Recently, deposits from the high terraces of the Rhine between Bonn (Germany) and Venlo (the Netherlands) have been dated to 750 ± 250 and 740 ± 210 ka using in situ-produced cosmogenic radionuclides (Dehnert et al., 2011).Similarly, Meuse terrace deposits in the Liège area (Romont, Belgium) that are generally assumed to correspond with the series of main terrace deposits, have been dated to 725 ± 120 ka using the same technique (Rixhon et al., 2011).
During the Middle Pleistocene, the hydrography of northern Belgium drastically changed due to the "opening" of the English Channel (Vandenberghe and De Smedt, 1979;Fig. 2).Various studies (Gibbard, 2007;Gupta et al., 2007;Toucanne et al., 2009) link the opening of the English Channel to the catastrophic drainage of a large proglacial lake during marine isotope stage 12 (MIS 12), approximately 450 ka ago (Elsterian).The 450 ka event triggered the formation of a buried paleo-channel system known as the Flemish Valley, with extensions towards the south and the east (Tavernier and De Moor, 1974).The Nete catchment is generally considered to be the eastern extension of the Flemish Valley.At present, the Campine Plateau is a landform that markedly stands out with respect to its surroundings.It is a fluvial terrace covered by coarse gravelly Meuse deposits in the south and southeast and sandy Rhine deposits in the north (Fig. 3).The sediments have proven to feature a periglacial paleoenvironment and were deposited by braided rivers (Paulissen, 1973 and1983).The Campine Plateau can be considered a classical case of relief inversion, given its prominent position in the landscape (Paulissen, 1983;Fig. 4).Undoubtedly, the area west of the Campine Plateau experienced prolonged phases of erosion and denudation after the Rhine had left the region, around 0.5 Ma at the latest (Fig. 4b; Beerten et al., 2017).

Cosmogenic radionuclide profiling
Cosmogenic radionuclides (CRN) allow us to quantify geomorphological processes over time spans covering the last 2 Myr.In this study, we use the concentration vs. depth profile of a single in situ-produced CRN ( 10 Be) to constrain the post-depositional denudation rate, E (m Myr −1 ) of the fluvial terrace.The accumulation of CRN, N total (z, t) (atoms g −1 ), in an eroding surface can be described by a mathematical function composed of two terms that represent the inherited CRN concentration of the fluvial sediment, N inh (atoms g −1 ), and the post-depositional production of CRN, N exp (z): where E is expressed in cm yr −1 (m Myr −1 × 10 −4 ), t (years) is the exposure age, λ (1 yr −1 ) the decay constant (λ = ln 2/t 1/2 ), z 0 the initial shielding depth (z 0 = E × t), ρ (g cm −3 ) the density of the overlying material and i (g cm −3 ) the attenuation length.The production rate of CRN, P i (z) (atoms g −1 yr −1 ), is a function of the depth, z (cm), below the surface: The subscript "i" indicates the different production pathways of in situ-produced 10 Be via spallation, muon capture and fast muons following Dunai (2010).In this study, the relative spallogenic and muogenic production rates are based on the empirical muogenic-to-spallogenic production ratios established by Braucher et al. (2011), using a fast muon relative production rate at SLHL of 0.87 % and slow muon relative production rate at SLHL of 0.27 %.The effective attenuation length is here equal to the apparent attenuation length as the depth profile was taken on a horizontal surface.The effective attenuation length for the sampling position was obtained using Table 4 in Marrero et al. (2016), and equals 152 g cm −2 .For fast and stopped muons, the at-Table 1. Analytical results from the in situ-produced 10 Be analysis.The depth profile is located at 50.95 • N and 5.63 • W at an altitude of 45 m.A sea-level, high-latitude (SLHL) production rate of 4.25 ± 0.18 atoms g −1 yr −1 was used, which represents the regionally averaged SLHL production for Europe (Martin et al., 2017).Refer to the main text for more information on the methodology used.tenuation length was set at resp.1500 and 4320 g cm −2 following Braucher et al. (2011).Production rates were scaled following Stone (2000) with a sea level high-latitude production rate of 4.25 ± 0.18 atoms g −1 yr −1 (Martin et al., 2017).The latter represents the regionally averaged SLHL production rate for Europe.The bulk density, ρ, of the studied fluvial sediment was set to 1.7 g cm −3 , which is consistent with the average value of upper Neogene and Quaternary sediments in the region (Beerten et al., 2010).A half-life of 1.387 ± 0.012 × 10 6 yr was used for 10 Be following Cmeleff et al. (2010).The CosmoCalc add-in for Excel was used to calculate the scaling factors.Given the flat topography of the Campine Plateau, topographic shielding was negligible and therefore not corrected for (Norton and Vanacker, 2009).

Sampling and analytical methods
The depth profile was sampled in a sand pit (SCR-Sibelco NV) on the northwestern edge of the Campine Plateau (Fig. 4a and b).The altitude of the sampling spot is ca.47 m (Tweede Algemene Waterpassing), while the crest of the plateau further east reaches an altitude of ca.48 m.The almost 4 m thick sequence is composed of medium-grained quartz-rich fluvial sands, overlain by a thin layer (35 cm thick) of fine-grained aeolian sand (Fig. 5).Detailed grainsize characteristics of the fluvial sand are given in Fig. 6.Note that sample depth is given with reference to the top of the fluvial sands.The lowermost unit A consists of medium sand with mode and median in the range between 250 and 500 µm, while a significant portion of grains coarser than 500 µm is present.Unit B is finer with a median grain size of ca.250 µm and virtually no coarse sand (i.e., > 500 µm).Unit C consists of coarse sand (median grain size more than 500 µm) with a significant amount of fine gravel fragments.The next unit (E) is the finest unit of the sequence, with mode From the depth profile, 10 samples were collected for CRN analysis at depths ranging from 45 to 355 cm below the surface, from which 9 were analyzed.Samples were more or Figure 6. 10 Be concentration profile and results of the grain size analysis.Note that the elevated 10 Be concentrations belong to samples that were analyzed using a smaller grain size fraction than the other samples (i.e., 250-500 µm instead of 500-1000 µm).These are indicated by pale gray dots.Depth is given relative to the uppermost sample.less evenly spread out over the sequence, although the sampling density was higher towards the top (Table 1).Samples were taken as bulk samples of 1.5 kg, over a depth interval of 10 cm.Samples were sieved, and the 500-1000 µm grain size fraction was used for sample preparation, except for the finegrained sand samples MHR-II-04 and MHR-II-06, where the 250-500 µm fraction had to be used.
Samples were prepared at the University of Louvain Cosmogenic Isotope Laboratory (Louvain-la-Neuve).In situproduced 10 Be was extracted from purified quartz using standard separation methods described in von Blanckenburg et al. (1996) and Vanacker et al. (2007b).Two blanks were processed with the nine samples.Approximately 200 µg of 9 Be carrier was added to blanks and samples containing 30 to 35 g pure quartz.The 10 Be / 9 Be ratios were measured in BeO targets with accelerator mass spectrometry on the 0.6 MV Tandy at ETH Zurich (Kubik and Christl, 2010).The ratios were normalized to the ETH secondary in-house standard S2007N with a nominal value of 10 Be / 9 Be of 28.1 × 10 −12 (Kubik and Christl, 2010), which is in agreement with a half-life of 1.387 Myr (Chmeleff et al., 2010).Samples are corrected for the number of 10 Be atoms in their associated blanks.The analytical uncertainties on the 10 Be / 9 Be ratios of sample and blank are then propagated into the 1σ analytical uncertainty for nuclide concentrations.

Inverse problem
To acknowledge that measurements and modeling errors are inevitable, the inverse problem is commonly represented by the stochastic relationship given by where d = (1, . .., N ) ∈ R N , N ≥ 1 is the measurement data, F is a deterministic forward model with parameters, x, and the noise term, e, lumps measurement and model errors.
Inversions were performed within a Bayesian framework, which treats the unknown model parameters x as random variables with posterior probability density function (PDF), p (x|d), given by where p (x) denotes the prior distribution of x and L (x|d) ≡ p (d|x) signifies the likelihood function of x.The normalization factor p (d) = p (d|x) p (x) dx is obtained from numerical integration over the parameter space so that p (x|d) scales to unity.The quantity p (d) is not required for parameter inference.Unless stated otherwise, in the remainder of this study we will focus on the unnormalized posterior p (x|d) ∝ L (x|d) p (x).
If we assume the residual errors, e, to be normally distributed, uncorrelated and with unknown constant variance, σ 2 e , the log-likelihood function can be written as For numerical stability, it is, however, often preferable to work with the log-likelihood function, l (x|d), instead of L (x|d): Earth Surf.Dynam., 5, 331-345, 2017 The variance of the residual errors, σ 2 e , can be fixed beforehand or sampled jointly with the other model parameters x.Note that by fixing σ 2 e to a known measurement error, σ 2 m , one implicitly assumes that the model is able to describe the observed system up to the observation errors.This might not be realistic in environmental modeling, where models are always fairly simplified descriptions of a more complex reality.In this work, we therefore jointly infer σ e with x.This accounts for both measurement and model errors, under the assumption that both types of errors obey a zero-mean uncorrelated and homoscedastic normal distribution.

Markov chain Monte Carlo sampling
The inference seeks to estimate the posterior parameter distribution of the model parameters, p (x|d).As an exact analytical solution for p (x|d) is not available, we resort to Markov chain Monte Carlo (MCMC) simulation to generate samples from this distribution.The basis of this technique is a Markov chain that generates a random walk through the search space and iteratively finds parameter sets with stable frequencies stemming from the posterior PDF of the model parameters (see, e.g., Robert and Casella, 2004, for a comprehensive overview of MCMC simulation).In practice, the MCMC sampling efficiency strongly depends on the assumed proposal distribution used to generate transitions in the Markov chain.In this work, the state-of-theart DREAM (ZS) (ter Braak and Vrugt, 2008;Vrugt et al., 2009;Laloy and Vrugt, 2012) algorithm is used to generate posterior samples.A detailed description of this sampling scheme including convergence proof can be found in the cited literature and is thus not reproduced here.Note that the considered CRN data inversion is a fairly simple problem (the model in Eqs.1-2 is quick and well-behaved, whereas both the parameter and measurement data spaces are rather low-dimensional).The use of DREAM (ZS) will become even more attractive when considering larger parameter dimensionality.

Prior distribution
The prior PDF is a key element of Bayesian inference.This distribution encodes the available prior information about the parameters and balances the effect of the likelihood function on the posterior PDF (Eq.4).We assumed the individual prior parameter PDFs to be independent: with N p = 4 the dimension of x.
Based on the current geological knowledge of the region, we specified a truncated Gaussian prior distribution for E, with mean of 30 m Myr −1 , standard deviation of 30 m Myr −1 and range of [0, 60] m Myr −1 .The upper bound of 60 m Myr −1 is based on (1) geomorphological evidence presented in Fig. 4b, using the altitude difference between the Campine Plateau and the adjacent Kleine Nete floodplain, and (2) the youngest possible age for the Rhine sediments.The lower bound of 0 m Myr −1 corresponds to the scenario where the Campine Plateau is a residual relief due to erosion resistance of the covering sediments.Overall, the resulting p (E) is sufficiently vague to avoid over-constraining the inversion while nevertheless discouraging the search to pick up values close to the boundaries that are considered to be (a priori) less likely.A uniform prior distribution in the range [0, 1] Myr was selected for t.This is based on the presumed burial age of the Rhine sands covering the Campine Plateau (between 0.5 and 1 Ma) together with geologic evidence on the evolution of the Scheldt Basin and in particular the Nete catchment after the opening of the English Channel (0.45 Ma, see Sect. 2).In addition, we put a uniform prior PDF with range [1, 35] m on the product E × t.This limits the total erosion that can possibly be inferred from the measurement data to 35 m, as we expect that 35 m of total erosion on top of the Campine Plateau is an absolute maximum.When adding the thickness of Rhine sediment that covers the top of the Campine Plateau (i.e., 5 to 10 m) to the maximum total erosion, we obtain a maximum initial thickness of 40 to 50 m, which corresponds to the thickness of Rhine deposits that are preserved in the deepest part of the Roer Valley graben (Beerten, 2006;Deckers et al., 2014).The minimum amount of total erosion is set to 1 m, given the altitude of the sampling position which is slightly (ca. 1 m) lower than the crest of the plateau.The prior distribution for the inherited 10 Be concentration, N inh , was assumed to be uniform.This is because solutions for N(z, t) in Eq. (1) always converge to N inh as z tends to infinity, while no N inh measurements are available at z larger than 3.5 m.The N inh parameter was therefore allowed to vary uniformly between 1 × 10 4 atoms g −1 and 9 × 10 4 atoms g −1 .The upper bound of 9 × 10 4 atoms g −1 is consistent with the highest inherited concentration measured in the profile (Table 2; the inherited concentration cannot be higher than this value).The lower bound was somewhat arbitrarily set at 1×10 4 atoms g −1 , given the fact that zero inheritance is considered to be very unlikely.Lastly, a so-called Jeffreys (1946) prior of the form p (σ e ) ∝ 1/σ e was used for σ e .This classical choice basically means that one wants to achieve σ e values that are both as small as possible and large enough to be consistent with the data misfit.

Comparison with the Bayesian approach in CRONUScalc
For comparison, we implemented the approach taken by Marrero et al. (2016) in CRONUScalc (CR).For brevity, in the reminder of this paper we will refer to this CRONUScalc Bayes method as CRB.Similarly as with our approach, CRB seeks to estimate p (x|d) using Eq. ( 4).However, CRB does not generate a set of samples with frequencies stemming www.earth-surf-dynam.net/5/331/2017/Earth Surf.Dynam., 5, 331-345, 2017 The need to evaluate p (d) using trapezoidal integration implies that CRB cannot use l (x|d) but requires using L (x|d) instead.The L (x|d) formulation taken by CRB is similar to Eq. ( 5), except for two aspects.First CRB includes the more general heteroscedastic case as well.In other words, the residual errors, e = [e 1 , . .., e N ], can have different variances, σ 2 e 1 , . .., σ 2 e N .Second and most important, CRB sets the σ e i to analytical measurement errors, σ m i .Using similar notations as in Marrero et al. (2016), CRB uses where χ 2 is a weighted sum of squared residuals (WSSR) also referred to as the chi-square statistic.In this study σ m is assumed to be constant: σ m 1 , . .., σ m N = σ m Equation ( 10) therefore reduces to Eq. ( 6) but with the standard deviation of the residual errors, σ e , fixed to σ m .Thus CRB considers three parameters: E, t and N inh .With respect to the associated 3-D prior distribution, p (x), we used the assumptions as for our approach (see Sect. 3.3.3).
As stated earlier, no matter whether ones uses a constant σ m or a different σ m i for each residual, e i , fixing the standard deviation(s) of the residual errors to the standard deviation(s) of the measurement errors implicitly assumes that the model can fit the concentration data within the standard deviation(s) of the measurement errors.If this assumption is not met, then the resulting p (x|d) estimation will be biased towards underestimating uncertainty.For the case of constant σ e and σ m in Eqs. ( 6) and ( 10), respectively, the solution that maximizes p (x|d), or maximum a posterior solution (MAP), should have a root mean square error, RMSE that is close to σ e (Eq.10) or σ m (Eq.6).Otherwise, the cho-sen likelihood model will not be consistent with the actual data.

Predictive uncertainty intervals
A 95 % uncertainty interval for the simulated CRN concentrations can be calculated by drawing parameter sets, x, from p (x|d) and removing the 2.5 % largest and lowest values from the associated set of F (x) responses.If all prior assumptions about the residual error distribution are met, then this 95 % predictive uncertainty interval should encompass 95 % of the measurement data.

CRN measurement data
In general, there is a clear decrease in 10 Be concentration with depth, except for two samples (MHR-II-04 and MHR-II-06) which contain higher CRN concentrations (Table 1 and Fig. 6).It is striking that the CRN concentrations are consistently higher for the two samples where the finer (250-500 µm) grain size fraction was analyzed.Grain size-dependent 10 Be concentrations can point to differences in geomorphological process rates in the regions of sediment provenance as suggested by Carretier et al. (2015).Alternatively, the negative relationship between in situ-produced CRN and grain size might also result from non-stationary sedimentation rates, where samples from the fine-grained layers accumulated CRN during the final stage of the sedimentation cycle prior to a phase of non-deposition and/or steady state.Apart from samples MHR-II-04 and MHR-II-06, the CRN concentrations decrease non-linearly with depth, from (1.5 ± 0.02) × 10 5 atoms g −1 at 45 cm to (9.0 ± 0.2) ×10 4 atoms g −1 at 355 cm (Table 1 and Fig. 6).Because they are not compatible with the other profile data, samples MHR-II-04 and MHR-II-06 were excluded from the inversion, thereby leading to a measurement data set of seven concentrations.These measured concentrations are associated with analytical measurement errors, σ m 1 , . .., σ m 7 , in the range of 6 × 10 3 -7 × 10 3 atoms g −1 .
If we consider the end-member where erosion rates of the Campine Plateau are very low (E ≈ 0 m Myr −1 ), as one could assume from its geomorphic setting as an inverted topography, the difference between N total (z) and N inh gives the N exp (z), the concentration of cosmogenic nuclides that is produced at depth z after deposition of the Rhine sands.The apparent exposure age of the surface, t, can then be reconstructed following Eq.(1).By doing so, we obtain an apparent exposure age of 21.5 ± 1.5 ka, which strongly contradicts the chronostratigraphical age estimates of the fluvial deposits that cover the Campine Plateau that range between 0.5 and 1 Ma (see Sects. 2 and 3.3.2).We advocate that postdepositional erosion has strongly altered the 10 Be signature of the upper layers of the Rhine sediments at the study site.

Posterior distribution derived using the approach proposed
We ran DREAM (ZS) for a (serial) total of 150 × 10 3 model evaluations.The marginal posterior distributions of the four sampled parameters (including the standard deviation of the residual errors, σ e ) are depicted in Fig. 7, while bivariate posterior scatter plots together with iso-density contour lines are presented in Fig. 8.The erosion rate, E, is relatively well resolved with a clear mode around 39 m Myr −1 (Figs.  2 and Fig. 8b) with a linear correlation coefficient of 0.67.The t posterior distribution shows a weakly expressed mode between ca.50 and 200 ka (Figs.7b, 8a and c).Also, t does not correlate with other sampled parameters, except for N inh (Table 2 and Fig. 8c) with which a linear correlation of −0.38 is observed (Table 2).The t parameter is therefore left largely unresolved by the inversion.
The N inh parameter shows a clear mode around 8.4 × 10 4 atoms g −1 (Figs.7c, 8b and c), which is lower than the lowest measured value in the profile of about 9.09 × 10 4 atoms g −1 (Table 1).As mentioned already, a moderately large dependence on E is observed (Fig. 8b and Table 2).The σ e parameter shows an approximately log-normal marginal posterior with a clear mode around 1 × 10 4 atoms g −1 (Fig. 7d).This is consistent with the RMSE values achieved between measured and simulated 10 Be.Indeed, the MAP solution induces an RMSE of 9.8×10 3 atoms g −1 .Notice that with values ranging between 6 × 10 3 and 7 × 10 3 atoms g −1 , the analytical measurement errors are 1.4 to 1.7 times smaller than the RMSE of the MAP solution.This illustrates the effect of model errors.If the model were perfect, the MAP solution would have been associated with an RMSE that is within the measurement error range of 6 × 10 3 − 7 × 10 3 atoms g −1 (Sect.4.1).Figure 9 presents the 95 % uncertainty interval associated with model predictions.This interval brackets five (71 % of the data) or six (86 % of the data) out of seven observations, depending on how the BE-MHR-II-03 measurement data point, which is located at the limit of the uncertainty band, is interpreted (Table 1).With only seven data points it is impossible to further assess the accuracy of the 95 % uncertainty band displayed in Fig. 9. Overall, it seems reasonably consistent.

Comparison against the approach taken in CRONUScalc
For CRB we sampled over an evenly spaced 3-D grid with upper and lower limits defined in Sect.3.3.3and 60 grid divisions in each dimension.This leads to a total of 216 × 10 3 model evaluations, which is a little bit more than what was used in our proposed approach (150×10 3 ).Moreover, we assumed a constant measurement error: σ m 1 , . .., σ m N = σ m = 7×10 3 atoms g −1 .Using a constant rather than variable measurement error is fully justified here because (1) the analytical error range is only between 6 and 7 × 10 3 atoms g −1 (Sect.4.1), and (2) as shown in Sect.4.2.1 the model cannot fit the data up to the maximum measurement error of 7 × 10 3 atoms g −1 anyway.The marginal posterior distributions of the three sampled parameters are presented in Fig. 10.The CRB finds similar modal or MAP values to our approach (compare Fig. 7ac with Fig. 10a-c).The t posterior distribution obtained by CRB is very close to that derived using our approach, except for a narrower peak around the MAP.For E and N inh , the posterior distribution obtained by CRB is a narrower version of that derived using our approach (compare Fig. 7a with 10a and Fig. 7c with 10c).This is caused by the use of σ m = 7 × 10 3 atoms g −1 in the likelihood function (Eq.10).The latter generally induces a more peaky likelihood (and consequently narrower posterior density) than our approach for which σ e values in the range shown in Fig. (8d) are sampled.Since, similarly as with our approach, the RMSE of the MAP solution derived using CRB is approximately 9.8 × 10 3 atoms g −1 , the CRB likelihood function is actually too narrow to be consistent with the achieved data misfit.Here, this leads to a moderate underestimation of uncertainty.This becomes more apparent in the resulting predictive uncertainty intervals (Fig. 11).Indeed the 95 % uncer-tainty band only brackets four of seven data points, that is, 57 % of the observations.Lastly, it is worth noting that the fact that the distributions presented in Fig. 7 are less smooth than those shown in Fig. 10 is due to the different natures of grid-based sampling and MCMC.A given bin in Fig. 7 is made of 3000 samples that are drawn from the posterior PDF by the MCMC, while each bin in Fig. 10 corresponds to a single (central-bin) point of the sampled lattice.

Accounting for model errors in CRONUScalc
A simple fix for making the CRB uncertainty estimates more consistent in the presence of model errors is as follows: I. Identify the minimum RMSE over the sampled lattice, plug it as an estimate of σ e in Eq. ( 6) and compute the posterior density of each lattice point.
II. Check whether the resulting MAP solution has an RMSE that is close to the fixed σ e .These two values will obviously be equal if uniform priors are used, but may not necessarily be similar otherwise.
III.If the previous condition (II) is satisfied, then proceed with the inference.Otherwise, set σ e to the RMSE of the MAP solution, recompute the posterior density of each lattice point and go back to II.
For the considered case study, this procedure expectedly leads to fixing σ e to about 9.8×10 3 atoms g −1 .This results in increased posterior parameter and predictive uncertainty that approach that derived using our approach in Figs.7 and 9.However, these uncertainty estimates remain slightly smaller than those displayed in Figs.7 and 9.This is because rather than fixing σ e , our approach infers its complete posterior distribution given the information content of the measurement data.

Discussion of the obtained erosion rate estimate
In Fig. 12, erosion rates for outcrops, as published by Portenga and Bierman (2011), are shown together with the mean erosion rate obtained in the present study.The global erosion data are based on surface samples (thickness ranging between 0.5 and 8 cm) from a variety of bedrock, including igneous, metamorphic and sedimentary rocks, and various climate-tectonic settings.Generally speaking, outcrop erosion rates from Portenga and Bierman (2011) tend to be lower than the 39 ± 8 m Myr −1 determined for the northwest- ern part of the Campine Plateau.This may be a reflection of the fact that the Portenga and Bierman (2011) dataset is partially based on bedrock samples which are more resistant to erosion than unconsolidated sediment.Nevertheless, in a western European context, the erosion rate that we report for the northwestern Campine Plateau seems to be fairly high for a fluvial terrace.For comparison, the Meuse younger main terrace (YMT) near Liège does not show any signs of postdepositional erosion following Rixhon et al. (2011).Probably, the coarse-grained and slightly consolidated nature of the Meuse gravels can be put forward as an explanation.Similarly, Dehnert et al. (2011) reported that the high terraces of the Rhine in Germany and the Netherlands were eroded by only 1-3 m, and that the loess cover presumably protected the Rhine sands from significant erosion soon after deposition.An alternative explanation for the relatively high erosion rate found in the current study may be the proximity of the North Sea, and low base level during glacial periods.In contrast to the Meuse and Rhine, the Scheldt Basin, to which the northwestern Campine Plateau belongs today, developed in response to the sudden base level lowering as a result of the opening of the English Channel, ca.450 ka (see Sect. 2).An important feature of the Scheldt Basin is the Flemish Valley, a buried river system.The sudden base level lowering may have caused a regressive erosion wave penetrating into the hinterland, shaping the Flemish Valley and its eastern extension, i.e., the Nete catchment, and cause increased erosion rates in this distal part of the Scheldt Basin.The posterior distribution for t, showing an increasing probability for t<0.5 Myr, and peaking between 200 and 50 ka supports this hypothesis.In the case of non-stationary erosion, it remains unclear to which extent the erosion rate can be used to infer the total amount of erosion at the site.In this study, we used Earth Surf.Dynam., 5, 331-345, 2017 www.earth-surf-dynam.net/5/331/2017/E × t as a joint prior for total erosion for which a lower and upper limit of 1 and 35 m was set, and our results show an erosion estimate distribution as given in Fig. 13.The posterior for E × t is poorly resolved, which is mainly caused by the poorly resolved posterior for t.
Our erosion estimate for the top of the northwestern Campine Plateau asks for a revision of regional landscape evolution models.Firstly, our results suggest that the total amount of Rhine sediment would have been larger than what can be observed today in the quarries (Fig. 4); this should be taken into account when correlating Rhine sediment from the Campine Plateau with that in the Roer Valley graben.Secondly, they indicate that post-depositional fault movement along (segments of) the Feldbiss fault as derived from the stratigraphy of Rhine sands should be considered to be a minimum (Fig. 3).Thirdly, the amount of post-depositional erosion west of the plateau (Fig. 2), as can be observed from present-day altitude differences (northwestern Campine Plateau vs. Nete Valley; Fig. 4b), should be regarded as a minimum erosion value.
In future work, we plan to consider more parameters within the Bayesian inversion when new and more densely sampled profiles become available.Increasing parameter dimensionality is straightforward with our MCMC sampling, but may quickly become intractable with the pure grid-based approach taken in the CRONUScalc program (Marrero et al., 2016).This limitation also holds for plain MC simulation as done by Hidy et al. (2010).Moreover, in an attempt to account for differences in geomorphological process rates in the regions of sediment provenance or non-stationary sedimentation rates, resulting in grain-size-dependent 10 Be con-centrations, our MCMC inversion could be combined with a distributed numerical forward modeling approach instead of the currently used analytical solution.

Conclusions
We inverted an in situ-produced 10 Be concentration depth profile within a Bayesian framework from the northwestern Campine Plateau (NE Belgium) to infer the average longterm erosion rate, surface exposure age and the inherited 10 Be concentration in the profile.Compared to the state of the art in probabilistic inversion of CRN profile data, our inversion approach has two new ingredients: it (1) uses Markov chain Monte Carlo (MCMC) sampling, and (2) accounts for (under certain conditions) the contribution of model errors to posterior parameter and predictive uncertainty.We compared our approach to that taken in the CRONUScalc program for the considered case study.Both approaches are found to produce similar maximum a posteriori (MAP) values.Nevertheless, the method implemented in CRONUScalc also moderately underestimates uncertainty.A simple fix for making these uncertainty estimates more consistent in the presence of model errors is therefore proposed.For the studied fluvial terrace of the Rhine which today belongs to the Nete catchment (Scheldt Basin), the derived MAP post-depositional erosion rate is ca.39 ± 8 m Myr −1 (1σ ).This is fairly high compared to published outcrop erosion rate data in the Meuse and Rhine catchment, and elsewhere in the world.We believe that the unconsolidated and gravel-poor nature of the studied Rhine sediment, the absence of a protecting cover (such as loess) and possibly also headward erosion in response to sudden base level lowering around 450 ka are possible explanations.Our future work will try to better resolve the erosion rate together with several other uncertain parameters from MCMC inversion of dense two-nuclide concentration depth profiles.

Figure 1 .
Figure 1.Location of the Campine region within Europe and the European sand belt.

Figure 2 .
Figure 2. Structural map of northwestern Europe showing the Roer Valley graben faults and the Brabant and Rhenohercynian (Ardennes) massifs superimposed on a digital terrain model (DTM; GTOPO30; data available from the U.S. Geological Survey), with indication of large rivers (http://www.eea.europa.eu/data-and-maps/data/wise-large-rivers-and-large-lakes) and location of the Scheldt Basin, (dashed line), the Nete catchment (dotted line) and the Flemish Valley (solid line).The general paleohydrography of the Meuse and Rhine between 0.5 and 1.0 Ma is shown in colored lines.Headward erosion as an explanation for the development of the Nete catchment is indicated with a yellow arrow.

Figure 4 .
Figure 4. (a) Detailed DTM of the study area (Digitaal Hoogtemodel Vlaanderen II, DTM, raster, 1 m), with indication of the sampling location (white arrow).Note the regularly shaped sand quarries south of profile line A-A' which appear as depressions on the DTM.(b) Topographic cross section according to the profile line (A-A') shown in (a).The sampling location is schematically shown as a gray rectangle.

Figure 5 .
Figure 5. Photograph of the sampled profile with indication of sampling points, field codes, lithological units (A-G) and approximate profile depth.

Figure 7 .
Figure 7. Marginal posterior distributions of the four parameters sampled using our approach: (a) erosion rate, (b) exposure age, (c) inherited 10 Be concentration and (d) standard deviation of the residual errors.The blue bars denote the posterior PDFs and the red lines signify the associated prior PDFs.

Figure 8 .
Figure8.Selected scatter plots together with iso-density contour lines for the sampled posterior parameter distribution by our approach.Black dots are posterior parameter sets and the density increases with the line color ranging from red (lower density) to yellow (higher density).

Figure 9 .
Figure9.95 % predictive uncertainty interval (gray area) and associated MAP prediction (black line) derived using our approach.The red crosses represent the seven measurement data points used in the inversion.Depth is given in absolute depth.

Figure 10 .
Figure 10.Marginal posterior distributions of the three sampled parameters using the Bayesian approach taken in the CRONUScalc program: (a) erosion rate, (b) exposure age, and (c) inherited 10 Be concentration.The blue bars denote the posterior PDFs and the red lines signify the associated prior PDFs.

Figure 11 .
Figure 11.95 % predictive uncertainty interval (gray area) and associated MAP prediction (black line) derived using the Bayesian approach taken in the CRONUScalc program.The red crosses represent the seven measurement data points used in the inversion.Depth is given in absolute depth.

Figure 12 .
Figure 12.Frequency distribution of outcrop erosion rates published by Portenga and Bierman (2011), with indication of the mean value obtained for the northwestern Campine Plateau (this study).

Figure 13 .
Figure 13.Derived posterior distribution for total erosion (denudation) at the study site (E × t).
Sample field code Relative depth Absolute depth Quartz Be carrier 10 Be/ 9 Be *The relative depth is given as depth below the uppermost sample.n/a: not applicable.

Table 2 .
Posterior linear correlation coefficients between the four parameters sampled using our approach.