Articles | Volume 4, issue 2
Geochronology, 4, 713–730, 2022
Geochronology, 4, 713–730, 2022
Research article
15 Dec 2022
Research article | 15 Dec 2022

Constraining the aggradation mode of Pleistocene river deposits based on cosmogenic radionuclide depth profiling and numerical modelling

Constraining the aggradation mode of Pleistocene river deposits based on cosmogenic radionuclide depth profiling and numerical modelling
Nathan Vandermaelen1, Koen Beerten2, François Clapuyt1,3, Marcus Christl4, and Veerle Vanacker1 Nathan Vandermaelen et al.
  • 1Georges Lemaître Centre for Earth and Climate Research, Earth and Life Institute, UCLouvain, Place Louis Pasteur 3, 1348 Louvain-la-Neuve, Belgium
  • 2Engineered and Geosystems Analysis, Waste and Disposal, Belgian Nuclear Research Centre SCK CEN, Boeretang 200, 2400 Mol, Belgium
  • 3Fonds de la Recherche Scientifique FRS-FNRS, Brussels, Belgium​​​​​​​
  • 4Laboratory of Ion Beam Physics, Department of Physics, ETH Zurich, Zurich, Switzerland

Correspondence: Nathan Vandermaelen ( and Veerle Vanacker (


Pleistocene braided-river deposits commonly represent long periods of non-deposition or erosion that are interrupted by rapid and short aggradation phases. When dating these sedimentary sequences with in situ-produced cosmic radionuclides (CRNs), simple concentration depth profiling approaches often fall short, as they assume that the alluvial sedimentary sequence has been deposited with a constant and rapid aggradation rate and been exposed to cosmic radiations afterwards. Numerical modelling of the evolution of CRNs in alluvial sequences permits one to account for aggradation, non-deposition and erosion phases and can simulate which scenarios of aggradation and preservation most likely represent the river dynamics. In this study, such a model was developed and applied to a Middle Pleistocene gravel sheet (Zutendaal gravels) exposed in NE Belgium. The model parameters were optimised to the observed 10Be and 26Al concentrations of 17 sediment samples taken over a depth interval of 7 m that constitutes the top of a gravel sheet up to 20 m thick. In the studied sedimentary sequence, (at least) three individual aggradation phases that were interrupted by non-deposition or erosion can be distinguished, each interruption lasting  40 kyr. The age for the onset of aggradation of the upper 7 m of the gravel sheet was further constrained to 654-62+218 ka. This age, within error limits, does not invalidate previous correlations of the gravel sheet with the Cromerian Glacial B and Marine Isotope Stage (MIS) 16. The deposition of the entire gravel sheet likely represents more than one climatic cycle and demonstrates the importance of accounting for the depositional modes of braided rivers when applying in situ cosmogenic radionuclide techniques.

1 Introduction

In situ-produced cosmogenic radionuclides (CRNs, e.g. 10Be and 26Al) are now widely used to infer erosion rates and exposure time of depositional landforms and allow one to better constrain the long-term landscape evolution of the Quaternary (e.g. Hancock et al., 1999; Schaller et al., 2001; Hidy et al., 2018). To constrain the post-depositional history of fluvial deposits, depth profiles are often used. They consist in measuring the CRN concentration over a depth interval of several metres below the surface. The CRN concentration in the upper 2 to 3 m decreases exponentially with depth, and the shape of the CRN depth profile informs on the average erosion rate and the post-depositional age. Below 3 m, the CRN concentration asymptotically decreases to a value that is assumed to represent the CRNs conserved from previous exposure episodes, the pre-depositional inheritance (e.g. Siame et al., 2004; Braucher et al., 2009; Hidy et al., 2010).

Measured CRN concentrations can then be fitted to numerical model predictions via an optimisation process. A minimum of five samples from the same, undisturbed sedimentary sequence is often necessary to obtain reliable results for exposure and erosion rates (Braucher et al., 2009; Hidy et al., 2010; Laloy et al., 2017). The depth profile technique assumes that the aggradation process is continuous and negligible in duration compared to the post-depositional exposure and that the inheritance is negligible or constant (Braucher et al., 2009; Laloy et al., 2017). Successful applications of CRNs to date Quaternary deposits include (glacio)fluvial terraces (e.g. Rixhon et al., 2011, 2014; Xu et al., 2019), alluvial fans (e.g. Rodés et al., 2011) or glacial moraines (e.g. Schaller et al., 2009) that underwent negligible or constant erosion rates over the time of exposure (Braucher et al., 2009).

Depositional sequences can show discontinuous aggradation modes, limiting the applicability of classical CRN depth profiling. Examples exist of Pleistocene river deposits that consist of several sedimentary cycles (Mol et al., 2000; Vandenberghe, 2001; Lauer et al., 2010, 2020; Vandermaelen et al., 2022a). Between the aggradation of each sequence lies a potential phase of landscape stability or erosion, hereafter referred to as a hiatus. While short hiatuses in the deposition process are, in principle, undetectable from a depth profile, recent work by e.g. Vandermaelen et al. (2022a) showed that > 5000 years-long hiatuses leave a clear imprint on the CRN depth profile. Further development of the classical depth profile technique is necessary to account for multiple aggradation phases and modes when constraining the history of depositional landforms like braided-river deposits (Nichols et al., 2002, 2005; Balco et al., 2005; Dehnert et al., 2011; Rixhon et al., 2011, 2014; Rizza et al., 2019).

In this study, we evaluate whether it is possible to reconstruct the aggradational mode of Pleistocene braided-river deposits based on in situ-produced CRN data collected over a  7 m thick sedimentary sequence. We developed a numerical model that simulates the accumulation of cosmogenic radionuclides, 10Be and 26Al, in a sedimentary sequence and that accounts for deposition and erosion phases and post-depositional exposure. The model is applied to Pleistocene gravel deposits, the Zutendaal gravels. The 7 to 15 m thick gravel sheets of the Zutendaal Formation are found in NE Belgium and are assumed to be of Middle Pleistocene age (Paulissen, 1983; Gullentops et al., 2001; Beerten et al., 2018). The thickness of the deposits, availability of geochemical proxy data and good preservation make them an excellent candidate for this study on complex aggradation modes.

2 Material and methods

2.1 Accumulation of CRN over time

2.1.1 Principles of numerical model

The model simulates the buildup of (i) a sedimentary sequence, including phases of aggradation, non-deposition and erosion, and (ii) the in situ-produced cosmogenic radionuclide concentrations in the sedimentary column. The exposure time of the sequence corresponds to the time elapsed since the onset of the deposition (i.e. aggradation) of the oldest, bottommost layer. The model treats the exposure time as the sum of a discrete number of time periods of variable duration (kyr) (Fig. 1). During a time period, there is either deposition of sediments with a given thickness (cm) on top of the pre-existing column (aggradation phase), surface erosion of a given amount of sediments (cm) (erosion phase) or landscape stability (no addition nor removal of sediments). The model allows one to specify the lower and upper bounds of the duration of the time periods and the amount of aggradation or erosion. A time period characterised by either stability or erosion can be defined to represent not only an interruption in the aggradation process but also the (often much longer) post-depositional time. In the following sections, the total length of all aggradation, erosional or stable phases is referred to as the “total formation time (kyr)”, whereas the period of time after abandonment is cited as the “post-depositional time (kyr)”. The sum of the two thus represents the exposure time (kyr).

Figure 1Structure of the model. The total exposure time is divided in a number of time periods, during which the sediment column is building up (aggradation phase), eroding (erosion phase) or not changing (stability).


The functioning of the model is exemplified in Fig. 2, where an example is given for the buildup of a sedimentary column in three phases. During the first time period, T1, from t0 until t1, the deposition starts during an aggradation phase. Aggradation is then interrupted by an erosion phase of duration T2, lasting from t1 until t2. After erosion, a new aggradation phase of duration T3 occurs between t2 and t3. After that, the fluvial sequences are abandoned and preserved until now, i.e. the end of the exposure time. The depth variation of in situ-produced cosmogenic radionuclides in the sedimentary sequence shows the effect of the complex aggradation history (Fig. 2) with two superimposed CRN depth profiles. The lower CRN depth profile developed between t0 and t1 and was truncated during the erosion phase of T2. If the profile was buried at great depth (typically > 10 m) and shielded from cosmic rays, no further accumulation of CRN occurred after t2. The upper CRN depth profile developed since the onset of the T3 aggradation phase, and the buildup of CRN continued after abandonment of the sequence. Such a discontinuous aggradation mode creates a CRN concentration depth profile that cannot properly be explained by model descriptions of classical “simple” CRN depth profiles.

Figure 2Illustration of the effect of discontinuous aggradation on the depth profile of cosmogenic radionuclide concentrations. (a) The sedimentary sequence consists of two aggradation phases (T1 and T3) that are interrupted by erosion or stability during T2. (b) The in situ-produced CRN profile shows two superimposed classical CRN depth profiles. The lower part of the profile represents two depth profiles: the dashed line illustrates the CRN depth profile when T2 undergoes erosion (E> 0), and the solid line represents the CRN depth profile when T2 corresponds to a stability phase (E= 0).


2.1.2 Model equations

The production rate of in situ cosmogenic radionuclides at a given depth, z (cm), in a sedimentary deposit can be described as follows:

(1) P i z = P i z 0 e - z ρ Λ i .

Pi(z0) (in atoms per gram of quartz per year, later referred to as at. gqtz-1 yr−1) is the production rate of CRN (10Be or 26Al) at the surface, z=z0 (cm), via the production pathway i, denoting either spallation by neutrons or capture of fast or negative muons. The attenuation length, Λi (g cm−2), is a measure of the attenuation of CRN production with depth and was set to 160, 1500 and 4320 g cm−2 for the production by, respectively, neutrons, negative muons and fast muons (Braucher et al., 2011). The dry bulk density of material is written as ρ (g cm−3). The model predefines the sea level high-latitude (SLHL) production rate for 10Be at 4.25 ± 0.18 at. gqtz-1 yr−1 (Martin et al., 2017), and the value is then scaled based on latitude and altitude of the site, following Stone (2000). 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 % for 10Be and, respectively, 0.22 % and 2.46 % for 26Al.

The CRN concentration changes as function of time and depth, following Dunai (2010):

(2) C z , t = C inh e - λ t + i P i z λ + ρ ε Λ i e - ρ z 0 - ε t Λ i 1 - e - λ + ρ ε Λ i t ,

where Cinh (at. gqtz-1) is the concentration of inherited CRNs from previous exposure before or during transport to the final sink, λ (yr−1) is the nuclide decay constant, and ε is the erosion rate (cm yr−1). We used a half-life of 1387 kyr for 10Be (Korschinek et al., 2010; Chmeleff et al., 2010) and 705 kyr for 26Al (Nishiizumi, 2004).

The model simulates the CRN concentrations during the buildup of a sedimentary column and considers phases of aggradation, stability and erosion. The model is discretised in 1 cm depth slices. The aggradation or erosion rate (cm yr−1) is obtained by dividing the total thickness of the sediments deposited or eroded during one sedimentary phase by the duration of the phase. Then, the model calculates the corresponding thickness (cm) of the layer to be aggraded or removed per aggradation or erosion phase as a function of the aggradation or erosion rate. Aggradation phases are discretised in time steps of 1 kyr. The thickness of material to be aggraded is distributed equally over each time step of the aggradation phase. When the value is not discrete, the model keeps track of the remaining values and adds them to the thickness to be aggraded over the next time step. For every centimetre added on top of the column, the depth values are dynamically adjusted and change from z to z+1. Compared to earlier work by, e.g. Nichols et al. (2002) or Rizza et al. (2019), the advantage of our approach is the flexible setup of the model, whereby the user can tune the model complexity and its parameters easily as to adapt it to a specific study case.

At each time step, the concentrations of 10Be and 26Al along the depth profile are dynamically adjusted, taking into account the production or removal of CRNs during each phase (erosion or stability) or time step. The inherited and the in situ-produced CRN concentrations are corrected for the natural decay of 10Be and 26Al as a function of the remaining exposure time. The radioactive decay of CRNs becomes important for Middle to Late Pleistocene deposits and is explicitly accounted for, in contrast to earlier work by e.g. Nichols et al. (2002, 2005). The 26Al and 10Be in situ production rates and the inherited concentrations are predefined in the model. By default, the model assumes that all sediments arrive in their final sink with an inherited 26Al /10Be ratio equal to the surface production ratio which is set at 6.75 (Nishiizumi, 1989; Balco and Rovey, 2008; Margreth et al., 2016). The 26Al /10Be ratio of the inherited CRNs can be adjusted in the model to allow for simulations with a 26Al /10Be production ratio of 8.0 to 8.4 at depth, as reported by Margreth et al. (2016) and Knudsen et al. (2019). Such simulations could then represent the aggradation of material that is sourced by deep erosion by, e.g. (peri)glacial processes (Akçar et al., 2017; Claude et al., 2017).

2.1.3 Fitting model outputs to 10Be and 26Al observed data

We used the reduced χ2 value as the optimisation criterion. The null hypothesis stipulates that the probability of finding a χ2 value higher than the one that we measure is likely. In that case, the expected distribution cannot be rejected and is conserved as a possible solution. After each model simulation, the reduced χ2 value (Taylor, 1997) is derived following Eq. (3):

(3) χ 2 = 1 d i = 1 n Y i obs - Y i mod σ i 2 ,

where Yiobs-Yimod is the difference between observed and modelled 26Al and 10Be concentrations; σi is the standard error encompassing all process and analytical errors, and d corresponds to the degrees of freedom of the dataset that is equal to the number of observations less the number of unconstrained model parameters. For each simulation, the measured reduced χ2 and its associated probability of finding a reduced χ2 of higher value is reported. The null hypothesis was rejected at 0.05 significance level, and the parameters of the associated simulation are stored as possible solutions.

2.2 Study area

We selected a study site in the Zutendaal gravels, a fluvial deposit covering the western part of the Campine Plateau (Fig. 3a). The Campine Plateau is a relict surface standing out of the otherwise-flat Campine area, a characteristic lowland of the European Sand Belt. This landscape is primarily a result of periglacial, fluvial and aeolian processes related to glacial–interglacial climatic cycles that took place during the Pleistocene (Vandenberghe, 1995). It is bordered in the east by the terrace staircase of the Meuse Valley and in the northeast by tectonic features, the Feldbiss Fault Zone and the Roer Valley Graben (Fig. 3b). In the southwest, the Campine Plateau is bordered by a cryopediment shaping the transition to the Scheldt basin (Beerten et al., 2018, and references therein).

Figure 3(a) Location of the Campine Plateau (black dashed line) in the low-lying region of the Campine area (CA, delineated by an orange dashed line) on a DTM (digital terrain model; U.S. Geological Survey's EROS Data Center: GTOPO30, 1996), with indications of main faults (Deckers et al., 2014). The whole area belongs to the sandy lowlands of the European Sand Belt. The Campine Plateau stands out of its environment by > 50 m. (b) The Zutendaal gravels (delineated with a brown line) cover the southeastern region and the most elevated parts of the plateau (DTM from GDI-Vlaanderen, 2014). They were deposited by the Meuse River by the time it was flowing westward (red arrow, P1). Later on, this course was abandoned as the Meuse River moved eastward to form the present-day Meuse Valley (red arrow, P2).

The Zutendaal gravels were deposited by the Meuse River during the course of the Early and/or Middle Pleistocene (Beerten et al., 2018; Fig. 3b). By this time, the region corresponded to a wide and shallow river valley occupied by braided-river channels (De Brue et al., 2015; Beerten et al., 2018, and references therein). The Zutendaal gravels are structured as superposed units that possibly represent different aggradation phases related to various deposition modes of braided rivers (Paulissen, 1983; Vandermaelen et al., 2022a). Architectural elements that support the existence of individual aggradation phases include gravel bars and bedforms, channels, sediment gravity flows and overbank fines (Dehaen, 2021). Such an assemblage approaches the structure of shallow gravel-bed braided rivers, defined as Scott-type fluvial deposits by Miall (1996), but has an unusual presence of clay plugs. After deposition of the gravels, the Campine area was subject to erosion (Beerten et al., 2013; Laloy et al., 2017). The erosion-resistant cap of gravel deposits played an important role in the Quaternary landscape evolution of the Campine area: the Zutendaal gravels are now observed at the highest topographic position in the Campine landscape as result of relief inversion (Beerten et al., 2018), whereby the gravel sheet has been covered by wind-dominated Weichselian sands of the Opgrimbie Member, Gent Formation (Beerten et al., 2017). Optically stimulated luminescence dating indicates an age range between ca. 23 and ca. 11 ka, covering the Late Pleniglacial and Late Glacial (MIS 2, Marine Isotope Stage; Derese et al., 2009; Vandenberghe et al., 2009).

The age control on aggradation and post-depositional erosion is poor, but the onset of aggradation is commonly assumed not to be older than 1000 ka (van den Berg, 1996; Van Balen et al., 2000; Gullentops et al., 2001; Westerhoff et al., 2008). The onset of post-depositional erosion remains unknown, but post-depositional erosion did not start before 500 ka (van den Berg, 1996; Van Balen et al., 2000; Westerhoff et al., 2008). The specific duration and mode of aggradation of the Zutendaal gravels remain currently unresolved.

2.3 Sampling and laboratory treatment

We sampled an abandoned gravel pit at the geosite “Quarry Hermans” (Bats et al., 1995) in As (510029.10′′ N, 53546.19′′ E, Fig. 3), where the Zutendaal gravels are exposed over a thickness of about 7 m. At least a comparable thickness of the Zutendaal gravels is supposed to be present underneath the exposure. The gravel sheet is covered by 60 cm of coversands, whereby the top of the profile reaches an altitude of 85 m a.s.l. The section was described in the field, with annotations of grain size and sorting, sedimentary structures, and traces of chemical weathering, including oxidation (Vandermaelen et al., 2022a). These observations allowed the subdivision of the profile into six units (U1–U6, Fig. 4). Over the depth range of 70–660 cm, we took 37 bulk samples for grain size and bulk elemental analyses. Seventeen samples were processed for in situ-produced CRN analyses: 14 for 10Be to construct a depth profile and 3 for 26Al to analyse the 26Al /10Be ratio in the main sedimentological units of the profile. Given that the uncertainties on 26Al quantifications are larger than those for 10Be because of the necessity to perform additional measurements of stable 27Al on Inductively Coupled Plasma Atomic Emission Spectroscopy (ICP-AES), the number of 10Be samples outnumbers the 26Al samples.

Figure 4Field observations used to constrain simulations. The log illustrates six units, labelled from U1 to U6, that were defined based on grain size, contact geometry, sedimentary structure, sorting and weathering traces observed in the field. The maximum grain size of the units is shown on the top x axis (C: clay, Si: silt, VFS: very fine sand, FS: fine sand, MS: medium sand, CS: coarse sand, VCS: very coarse sand, VFG: very fine gravel, MG: medium gravel, CG: coarse gravel, VCG: very coarse gravel). They are covered by Weichselian, aeolian coversands (abbreviated UWS). U1, U3 and U4 present imbrication and represent crudely bedded gravels (Gh). U2, U5 and U6 represent horizontally bedded sands (Sh) or pebbly fine to very coarse sands (Sp). The granulometry is illustrated by the D50 (black line), and the limits of the shaded area define the D10 on the left and the D90 on the right of the black line. Panels (a)(d) represent illustrations of the different units. Panels (e)(f) represent very coarse gravels collected in U1 showing Mn coatings. Panel (g) illustrates weathered gravels from U3 showing Fe coatings.


Samples were processed for in situ cosmogenic 10Be and 26Al analyses, following Vanacker et al. (2007, 2015). Samples were washed, dried, and sieved, and the 500–1000 µm grain size fraction was used for further analyses. Chemical leaching with a low concentration of acids (HCl, HNO3 and HF) was applied to purify quartz in an overhead shaker. Later on, purified samples of 10–40 g of quartz were leached with 24 % HF for 1 h to remove meteoric 10Be. This was followed by spiking the sample with 9Be and by total decomposition in concentrated HF. About 0.200 mg of 9Be was added to samples and blanks. The beryllium in the solution was then extracted by ion exchange chromatography, as described in von Blanckenburg et al. (1996). Three laboratory blanks were processed.

The 10Be /9Be and 26Al /27Al ratios were measured using accelerator mass spectrometry on the 500 kV Tandy facility at ETH Zurich (Christl et al., 2013). The 10Be /9Be ratios were normalised with the in-house standard S2007N and corrected with the average 10Be /9Be ratio of three blanks of (1.3 ± 0.6) × 10−14. The analytical uncertainties on the 10Be /9Be ratios of blanks and samples were then propagated into the 1 standard deviation analytical uncertainty for 10Be concentrations. We then plotted the 10Be concentrations and their respective uncertainties as functions of depth below the surface.

We measured the 27Al concentrations naturally present in the purified quartz by inductively coupled plasma-atomic emission spectroscopy (Thermo Scientific iCAP 6000 Series) at the MOCA platform of UCLouvain in Louvain-la-Neuve, Belgium. The 26Al /27Al results measured at ETH Zurich were calibrated with the nominal 26Al /27Al ratio of the internal standard ZAL02, equal to (46.4 ± 0.1) × 10−12. We subtracted from each measurement a 26Al /27Al background ratio of (7.1 ± 1.7) × 10−15 (Lachner et al., 2014). We reported analytical uncertainties with 1 standard deviation, including the propagated error of 24 % coming from the background, the uncertainty associated with AMS counting statistics, the AMS external error of 0.5 % and the 5 % stable 27Al measurement (ICP-AES) uncertainty. The accuracy of the element chemistry was tested with reference material BHVO-2, and the analytical uncertainty was evaluated at < 3 % for major element concentrations and < 6 % for trace element concentrations (Schoonejans et al., 2016).

2.4 Scenarios to constrain the geomorphic history of fluvial deposits

We implemented four scenarios in our model. Each scenario consists of a succession of “n” time periods, referred to as Ti, that correspond to the interval [ti−1, ti], with i representing the limits of a time period in the geomorphic history. Parameters are listed in Table 1. The periods are characterised by a single geomorphic setting (i.e. E is erosion, A is aggradation, S is stability), whereby a given thickness of sediment is removed or aggraded (cm) over a certain time interval (kyr). Durations and thicknesses are sampled from a uniform distribution, whose upper and lower bounds are stated between square brackets in Table 1. Inheritance parameters were sampled from a normal distribution (with mean and standard deviation reported in Table 1) that is centred on the inheritance that was reported in previous studies on Quaternary Meuse deposits (Rixhon et al., 2011; Laloy et al., 2017). A uniform bulk density of 2.1 g cm−3 was used, based on bulk density measurements of 17 samples. The four scenarios are summarised in Fig. 5.

Table 1Description of four scenarios that are used to constrain the geomorphic history of fluvial deposits. The different periods are presented in the headers. Each period is characterised by aggradation (A), stability (S) or erosion (E); a specific duration; a thickness to aggrade or remove; and an inherent concentration. When uniform distributions are used, the lower and upper bounds of the interval are given between square brackets. When values are depicted in a normal distribution, the mean and standard deviation are given between round brackets.

Download Print Version | Download XLSX

Figure 5The four scenarios that are developed to represent the sedimentary sequence of the Zutendaal gravels in As. Grey arrows represent aggradation (upward arrow), erosion (downward arrow) or stability (horizontal arrow). The material that is removed by erosion is coloured in grey.


We ran each scenario 107 times with 6.75 as the inherited 26Al /10Be ratio and then again 107 times with a ratio of 7.40. The latter is a mean value between the 26Al /10Be ratio for production at the surface and the ratio observed at depth by e.g. Margreth et al. (2016). By varying the inherited 26Al /10Be ratio, we aim to account for a potential mix of sediments sourced by deep and surface erosion. A plot of kernel-density estimates using Gaussian kernels was then generated from the simulations that were considered to be possible solutions based on their associated reduced χ2 value. For each parameter, the value representing the highest density of solutions is given as the optimal model outcome, and the uncertainties are reported with 95 % confidence intervals (2σ).

Based on prior information on the evolution of the Zutendaal gravels, four constraints were defined:

  1. “Not longer than”. The total geomorphic history takes place over the last 1000 kyr (Beerten et al., 2018, and references therein). In consequence, any scenario for which the total duration exceeded 1000 kyr was automatically discarded.

  2. “Not shorter than”. In the area where the Zutendaal gravels are currently outcropping, deposition ended by 500 ka at the latest (Westerhoff et al., 2008). The abandonment of the terrace must be older than 500 ka.

  3. “Final thickness of every unit should correspond to the measured thickness”. After an aggradation phase, the thickness of the unit can decrease by erosion over the following time period until it matches the observed, present-day thickness.

  4. The two last time periods of any scenario should include 60 to 200 cm aggradation of Weichselian coversands (unit Weichselian coversands, abbreviated “UWS”), followed by an erosion phase until the thickness of the coversands matches the present-day thickness.

The model was run for each scenario using the constraints and parameter distributions of Table 1. Parameter values were attributed randomly, following Hidy et al. (2010). The two first scenarios, scenarios 1 and 2, represent the classical depth profile (Braucher et al., 2009): scenario 1 represents a long and constant post-depositional erosion phase, whereas scenario 2 represents a long stable phase that is followed by a recent, pre-Weichselian episode of rapid erosion. In scenario 1, the fluvial gravel sheet starts accumulating CRNs over a [500, 1000] kyr period that is characterised by constant erosion of [0, 500] cm. In contrast, in scenario 2, the fluvial sheet remains stable over T1 and then undergoes a rapid erosional phase of [0, 500] cm over 1 or 10 kyr (T2). The maximum erosion rate is thus 500 cm in 1 kyr, or 5 mm yr−1, and is based on the upper range of long-term erosion rates reported in literature (Portenga and Bierman, 2011; Covault et al., 2013). After erosion, the fluvial sheet is covered by Weichselian coversands (UWS) during a 1–10 kyr period. The sand deposit is then eroded until it reaches the present-day thickness of 60 cm. In both scenarios, the onset of the in situ CRN accumulation is concomitant with the beginning of the post-depositional period (Fig. 5).

In contrast to the instantaneous aggradation mode of scenarios 1 and 2, scenarios 3 and 4 consider a stepped aggradation mode and are based on ancillary data from geochemical proxies (Vandermaelen et al., 2022a). As illustrated in Fig. 5, the simulations start with unit U1–U2 in place at the bottom of the sedimentary sequence. This is followed by different phases of aggradation and erosion. The duration of any erosion phase is set to a maximum of 60 kyr so that two erosion phases account for a maximum of 120 kyr. This corresponds to the duration of a full glacial cycle, i.e. about 110 to 120 kyr, based on Busschers et al. (2007). The first phase of erosion of [0, 500] cm is thus simulated over a T1 period of [0, 60] kyr. This is followed by the T2 period with [185, 685] cm aggradation of unit U3. The minimum aggradation corresponds to the present-day thickness of U3 and takes place during a single step, 1 kyr, or 10 kyr (Table 1). During the T3 period, the thickness of U3 that exceeds 185 cm is eroded over [0, 60] kyr. The next aggradation phase, T4, is characterised by [275, 775] cm aggradation of units U4, U5 and U6 over 1 or 10 kyr. The minimum aggradation corresponds to their present-day thickness. After this phase, the post-depositional evolution of the fluvial sequence starts: in scenario 3, the thickness of U4, U5 and U6 that exceeds 275 cm is then slowly and constantly eroded over [500, 1000] kyr, corresponding to the T5 period. In scenario 4, there is a long phase of stability during the T5 period that is then followed by a phase of rapid erosion during the T6 period, lasting [1, 10] kyr. In both scenarios 3 and 4, the last two periods are similar to the aggradation and erosion of the Weichselian coversands specified in scenarios 1 and 2.

3 Results

3.1 In situ-produced CRN concentrations along the depth profile

The 10Be concentrations vary from 120 × 103 to more than 200 × 103 at. gqtz-1 (Table 2). The total uncertainties on the measured 10Be concentrations are below 7 %, with the exception of the lowermost sample (Heras-02). The observed CRN depth variation deviates from a simple exponential decrease of 10Be concentration with depth and points to a complex deposition history. The upper eight values (70 to 300 cm, corresponding to U6, U5 and U4) show an exponential decrease from (175 ± 8) × 103 to (122 ± 7) × 103 at. gqtz-1 (Fig. 6). There is an abrupt change in the concentration at 370 cm depth, corresponding to (135 ± 7) × 103 at. gqtz-1 measured at the top of the U3 unit, about 12 % higher than the sample taken at 300 cm depth. The following two samples in U3 show a steady decrease in 10Be concentration with depth. At the bottom of the profile, in the upper part of U1, a third local maximum of 10Be concentration is found. The value of (202 ± 8) × 103 at. gqtz-1 measured at 550 cm depth is the highest value that was measured in the profile and is 50 % higher than the average 10Be concentration measured in the overlying units. The samples taken in U1 show a steady decrease of 10Be concentration with depth, from (145 ± 8) × 103 at. gqtz-1 at 586 cm to (130 ± 30) × 103 at. gqtz-1 at 657 cm depth.

Table 2In situ-produced CRN concentrations. All reported uncertainties are 1 standard deviation (1σ).

Download Print Version | Download XLSX

Figure 6Panels (a) and (b) present the 26Al and 10Be concentration and (c) the 26Al /10Be ratio. The reported uncertainties are 1 standard deviation (1σ).


The three samples analysed for 26Al show a decrease of 26Al concentration with depth: from (934 ± 103) × 103 at. gqtz-1 at 197 cm depth (U5) to (734 ± 88) × 103 at. gqtz-1 at 370 cm depth (U3) and finally to (720 ± 122) × 103 at. gqtz-1 at 657 cm depth (U1). Considering the total uncertainties of 11 % to 17 % on the AMS and ICP-AES measurements, only the 26Al concentration of the uppermost sample is significantly higher than the two deeper ones (Fig. 6). Although based on a limited number of samples, the depth evolution of the 26Al concentrations differs from what is observed for the 10Be concentrations, as the 10Be concentration appears higher at 370 cm depth than at 197 cm depth.

The three 26Al /10Be ratios decrease with depth, with values of 7.41 ± 0.92 at 197 cm, 5.45 ± 0.73 at 370 cm and 5.54 ± 1.55 at 657 cm. The measured ratios are consistent with near-surface production (i.e. 26Al /10Be ratio of 6.75; Granger and Muzikar, 2001; Erlanger et al., 2012).

3.2 Optimal model fits

The optimal model fits for the scenarios representing the instantaneous aggradation mode (i.e. scenarios 1 and 2, Fig. 5) return a minimised reduced χ2 value above 11 and fail to represent the observed 10Be and 26Al data correctly. Best fits for these scenarios are also unsensitive to the inherited 26Al /10Be ratio. The optimal model fits for the scenario that consider a stepped aggradational mode and long and average erosion (i.e. scenario 3, Fig. 5) have a reduced χ2 value of 130 when using aggradation phases of 1 kyr and 147 using phases of 10 kyr. The goodness of fit does not improve when using an inherited 26Al /10Be ratio of 7.40 instead of 6.75 (Table 3).

Table 3Reduced χ2 values of the optimal model fits for scenario 1, 2, 3 and 4, with 26Al /10Be ratios of 6.75 and 7.40 and with 1 and 10 kyr durations of the aggradation phases. The star indicates whether the p value did not show a significant disagreement between observed and modelled CRN concentrations.

Download Print Version | Download XLSX

The scenarios that consider a stepped depositional history and a period of [500; 1000] kyr of landscape stability (scenario 4) show better optimal fits (Table 3). At the 0.05 significance level, the simulations can be accepted as possible solutions when the reduced χ2 value is below 1.83. With the inherited 26Al /10Be ratio of 6.75, the optimal simulations with aggradation phases of 1 and 10 kyr have a reduced χ2 value of, respectively, 1.55 and 1.44. The model fit improves when using inherited 26Al /10Be ratios of 7.40, with optimal reduced χ2 values of, respectively, 1.36 and 1.25. Given that the goodness of fit is similar for the simulations with aggradation phases of 1 or 10 kyr, we pooled all results in the kernel-density plots. Figure 6 resumes the results on the overall aggradation time of the entire sequence; the deposition age of the lowermost unit, U1; and the duration of hiatuses and corresponding surface erosion at the top of the U2, U3 and U6 layers. The optimal values are reported with their 95 % confidence intervals (2σ).

The hiatus at the top of the U2 unit is characterised by a duration of 44-7+15 kyr and an erosion amount of 17-17+56 cm in the simulations with inherited ratios of 6.75 (Fig. 7a) and 50-12+10 kyr and 23-23+56 cm in simulations with inherited ratios of 7.40 (Fig. 7b). The combination of hiatus duration and erosion amount results in an erosion rate of 4.60 and 3.86 m Myr−1 for simulations with inherited ratios of, respectively, 6.75 and 7.40. Although this erosion rate is 1 order of magnitude lower than the median global denudation rate (i.e. 54 m Myr−1) that was calculated from a set of 87 drainage basins (Portenga and Bierman, 2011), it is coherent with the lowest long-term incision rate reported for the Meuse catchment near Liège (Van Balen et al., 2000).

Figure 7Density plot of the possible model outcomes for scenario 4. Any simulation that returned a reduced χ2 smaller than 1.83 is considered to be significant and is included as a possible solution in the density plots. The density of significant solutions is shown by colour coding. The four upper panels (a, c, e, g) represent the parameters of simulations (n= 486) characterised by a 26Al /10Be of 6.75, and the lower panels (b, d, f, h) represent the parameters of simulations (n= 1695) with a 26Al /10Be of 7.40. For each parameter, the most likely value (i.e. the value with highest density of significant solutions) is reported along with its 95 % confidence interval (2σ). The 95 % confidence interval is derived from the 2.5 % and 97.5 % limits of the kernel cumulative density function.


For the hiatus at the top of U3, the highest density of solutions is observed around 42-35+18 and 45-42+14 kyr for simulations with an inherited ratio of, respectively, 6.75 and 7.40 (Fig. 7c, d), which is similar to the hiatus' duration established for the lower U2 unit. In contrast, the associated solutions for erosion are 1 order of magnitude higher than at the top of U1–U2, with values of 395-335+120 cm (94 m Myr−1) and 430-363+76 cm (96 m Myr−1) for inherited ratios of, respectively, 6.75 and 7.40. The hiatus on top of U6 encloses the time between the last aggradation of fluvial deposits and the Weichselian coversands. The highest density of possible solutions is found at 540-52+209 kyr for simulations with an inherited ratio of 6.75, and a similar value, i.e. 540-50+263 kyr, is found with a ratio of 7.40 (Fig. 7e). The erosion of the overburden that was present on top of U6 took place before the aggradation of Weichselian coversands and is estimated at 295-27+208 and 325-43+151 cm for inherited ratios of, respectively, 6.75 and 7.40 (Fig. 7f). Constant erosion of the overburden over the period of [500; 1000] kyr is not likely (i.e. scenario 3, Table 3), and a scenario with rapid erosion of the overburden results in significantly better model fit (i.e. scenario 4, Table 3). The current model setup fixed the time interval at [1, 10] kyr, which results in optimal values of 295 and 325 m Myr−1. Further research is needed to better constrain the timing and length of such an erosion phase, as they have a strong impact on the erosion estimates.

Based on the CRN age–depth modelling, the deposition age of the bottommost deposits of the exposure was constrained at 654-62+218 and 669-58+272 ka for inherited ratios of, respectively, 6.75 and 7.40 (Fig. 7g, h). The total formation time of the sedimentary sequence, including the aggradation phases and the sedimentary hiatuses, is estimated at 90-45+37 and 98-30+34 kyr for inherited ratios of, respectively, 6.75 and 7.40 (Fig. 7g, h). The top of the Zutendaal Formation that is exposed in As has a deposition age of 562-45+211 and 565-44+269 ka for inherited ratios of, respectively, 6.75 and 7.40.

4 Discussion

4.1 Aggradation mode of the Zutendaal gravels

The onset of aggradation of the Middle Pleistocene gravel deposits of the so-called Main Terrace of the Meuse in NE Belgium, the Zutendaal gravels, was commonly assumed to be situated between 500 and 1000 ka (van den Berg, 1996; Van Balen et al., 2000; Gullentops et al., 2001; Westerhoff et al., 2008). Here, by applying a numerical model to the CRN concentration depth profiles, the age for the onset of U1 aggradation is estimated to be 654-62+218 ka. This age, within error limits, agrees with previous correlations of the gravel sheet with the Cromerian Glacial B and Marine Isotope Stage (MIS) 16 (Gullentops et al., 2001). Around this time, the gravel sheet was already in the process of being built, and only the upper part of the Zutendaal gravels (potentially 15 to 20 m thick gravel sheet) is here exposed. Therefore, it cannot be excluded that the base of the (buried part of the) Zutendaal gravels corresponds with MIS 18.

The model simulations also confirmed the existence of a stepped aggradation mode, whereby phases of aggradation are alternating with phases of stability or erosion. At least three phases of aggradation with two intraformational hiatuses are identified: a first hiatus at the top of the U2 unit lasting 44-7+15 kyr and resulting in surface erosion of 17-17+56 cm; a second hiatus at the top of the U3 unit with similar duration, i.e. 42-35+18 kyr, but much higher surface erosion of 395-335+120 cm. The third and final aggradation phase (U4–U5–U6, ending at last at 562-45+211 ka) occurred before abandonment of the region by the Meuse River. The results are here reported for inherited ratios of 6.75 and aggradation phases of 1 and 10 kyr, but they do not differ significantly when using alternative inherited 26Al /10Be ratios (Table 3).

The total formation time of the sedimentary sequence exposed in As is estimated at 90-45+37 kyr. Given that only the upper 7 m are exposed in As (Gullentops et al., 2001), the deposition of the entire sequence of the Zutendaal gravels probably represents more than one climatic cycle. Slow aggradation leading to 10Be enrichment with depth (e.g. Nichols et al., 2002; Rixhon et al., 2014) is not directly observed in this sequence. However, it could eventually have occurred in the shallow and sandy unit, U2, that was not sampled for 10Be. Slow aggradation would further increase the total formation time. Such a prolonged formation time can occur in fluvial depositional systems in the absence of tectonic uplift and consecutive river downcutting (Sougnez and Vanacker, 2011). Only after the Meuse River abandoned its northwestern course (Fig. 1) did it develop a staircase of alluvial terraces in response to tectonic uplift of the Ardennes–Rhenish Massif (Beerten et al., 2018 and references therein).

According to simulations of scenario 4, an overburden remained in place from the abandonment time until it was removed by an erosion phase, which preceded the sedimentation of the Weichselian aeolian cover. Other scenarios could be tested, whereby the post-depositional history is characterised by alternating phases of aeolian sand deposition and subsequent removal of this loose material (Beerten et al., 2017). The aeolian sand deposits that cover the top of various Meuse River terraces in the northern part of the Meuse Valley are commonly assumed to be Weichselian in age (e.g. Vanneste et al., 2001; Derese et al., 2009; Vandenberghe et al., 2009) regardless of their morphological position and age (Paulissen, 1973). This being said, a few very thin Late Saalian aeolian sand deposits were found on top of the Saalian Eisden-Lanklaar terrace near the steep valley wall of the Meuse River east of the study site (Paulissen, 1973). This would imply that only the latest episode of sand movement is well preserved on the top of the terrace landscape and previous aeolian sand covers have been reworked during subsequent glacial–interglacial cycles that followed the deposition of the terrace deposits.

4.2 Added value of modelling complex aggradation modes

Studies using the classical CRN depth profile approach mainly envision delivering an exposure age of the surface, i.e. of the uppermost deposits only. By using numerical modelling, it becomes possible to reconstruct the aggradation time and mode of a sedimentary sequence based on the evolution of the CRN concentrations with depth. This case study in NE Belgium demonstrates that the total formation time of braided-river deposits, such as the Zutendaal gravels, may constitute nearly 20 % of the deposition age (Fig. 7g, h). This shows the importance of considering the aggradation mode of braided rivers when applying CRN techniques to > 3 m deep sedimentary sequences.

Firstly, this approach can account for the presence of hiatuses in the sedimentary sequence. Hiatuses in the aggradation process temporarily expose parts of fluvial sheets that would otherwise be buried at depth and partially shielded from CRN accumulation. This can create a positive offset in the CRN concentration depth profile, whereby the concentration at a given depth is higher than the true inheritance value (Fig. 6; Vandermaelen et al., 2022a). Such observations do not fit in the simple concentration depth distribution and are often classified as outliers or as results of differential inheritance (e.g. following the concepts reported by Le Dortz et al., 2011).

Secondly, by considering the aggradation mode, additional information on the sourcing of sediments can be extracted from the depth distribution of 26Al /10Be ratios. The 26Al /10Be ratios that are measured in fluvial sediments result from (i) the inherited 26Al /10Be ratios of the source material and (ii) the in situ CRN accumulation when the material is exposed to cosmic rays. Accounting for the changing depth of the sedimentary layers within the fluvial sheet may allow one to explain high 26Al /10Be ratios (i.e. > 6.75) measured at depth. For the case study in NE Belgium, the model fit improved slightly when an inherited 26Al /10Be ratio of 7.40 was used in the simulations (Fig. 7). Inherited 26Al /10Be ratios that are substantially higher than 6.75 often point to intense physical erosion in the headwater basins where material is sourced from deep erosion by e.g. (peri)glacial processes (Claude et al., 2017; Knudsen et al., 2019). The material that was buried several metres below the surface is then delivered to the fluvial system and breaks into smaller parts on its route to the final sink or during intermediate storage. Such deposits are then constituted of a mix of sediments characterised by different 26Al /10Be inherited ratios. In their final sink, the sediments will further accumulate CRN following the 26Al /10Be surface production ratio of  6.75 when exposed at (or close to) the surface, or they can maintain an 26Al /10Be isotope ratio above 6.75 when quickly buried to a depth above 300 g cm−2, where the relative production of 26Al to 10Be is larger due to muogenic production.

4.3 Trade-off between model complexity and sample collection

Optimisation methods like the reduced χ2 require that the number of observed data is larger than the number of free parameters (Hidy et al., 2010). There is thus a trade-off to make between the complexity of the model and the number of data points obtained by CRN analyses. In case of CRN concentration depth profiles, the number of samples is often limited by the capacity and financial constraints that are needed to process samples for CRN analyses. In the numerical model, the phases of aggradation are characterised by their duration, aggradation thickness and inherited 10Be concentration. They are followed by phases of stability or erosion with an unknown duration.

The inclusion of aggradation modes in CRN concentration depth profile modelling therefore requires more unconstrained parameters than the classical depth profile approach (Table 1). The minimum number of CRN measures that is required can be calculated as follows:

(4) N CRN > ( k × N aggrad ) + ( l × N erosion ) ,

where NCRN is the number of CRN measurements for individual samples; k and l are the number of unconstrained parameters for, respectively, each aggradation and erosion or stability phase. Further complexification of the model by e.g. including unconstrained parameters for sediment density or 26Al /10Be inherited ratio will further increase the required number of CRN observations. With a limited number of samples, the complexity of the model can be reduced by, e.g. keeping the CRN inheritance fixed. This is acceptable when the local minima of the CRN concentrations in a given sedimentary sequence are within 10 % of the mean of the lowest values. For braided-river deposits, it is also viable to fix the length of the aggradation phase to 1 or 10 kyr: most sedimentary sequences represent long periods of non-deposition or erosion interrupted by rapid and short-lived depositional events (Bristow and Best, 1993).

4.4 Chronostratigraphical implications

Based on the age of 654-62+218 ka obtained for the base of the sedimentary sequence in As, the Zutendaal gravels exposed at the geosite in As were most likely deposited during MIS 16, and the lower part maybe even during MIS 18, according to northwest European chronostratigraphical subdivision and correlation with the marine isotope record (Cohen and Gibbard, 2019). The uppermost part of the deposits (U6, Fig. 6) in As is dated at 562-45+211 ka, thereby yielding a youngest possible deposition age of 517 ka. This age would associate the uppermost deposits with MIS 14. However, it cannot be ruled out that, by that time, the Meuse had already shifted its course towards the east of the As sampling site (Van Balen et al., 2000). If this is true, the uppermost depositional sequence in As should correspond to material deposited by another local river rather than by the Meuse itself. It is possible that the Bosbeek occupied this part of the Campine Plateau after the Meuse had shifted to the east. A fossil valley floor of the Bosbeek River has been identified by Gullentops et al. (1993), only a few hundred metres away from the As sampling site.

The absolute dating of the Zutendaal gravels contributes to constraining the chronostratigraphical framework of the Campine Plateau and the Meuse River terraces. Several authors (e.g. Pannekoek, 1924; Paulissen, 1973) have formulated the hypothesis that the Zutendaal gravels may represent the Main Terrace in the part of the Meuse Valley downstream of Maastricht. However, the correlation between terrace fragments located at different locations along the Meuse River is still a subject of scientific debate. However, it remains interesting to compare the chronological framework of the Zutendaal gravels on top of the Campine Plateau with the younger Main Terrace levels of the Meuse River in the Liège area and the Ardennes. van den Berg (1996) used paleomagnetic techniques on the different sub-levels of the Sint-Pietersberg Terrace to obtain age estimates of the Main Terrace in the region of Maastricht. The age of the uppermost sub-level was estimated at 955 ka by van den Berg (1996) and at 720 ka by Van Balen et al. (2000) based on the same data. Both studies agree on an age of 650 ka for the next sub-level. Rixhon et al. (2011) dated the younger Main Terrace of the Meuse River in the locality of Romont,  25 km upstream of our sampling position, using in situ-produced 10Be depth profiles. They obtained an age of 725 ± 120 ka (mean ± 1 SD). Although this age is somewhat older than the top of the Zutendaal gravels that we dated with CRN depth modelling (i.e. 562-45+211 ka, optimal solution with 95 % CI), it is not inconsistent, as it overlaps with the CRN optimal solution considering the 95 % confidence interval.

5 Conclusion

The aggradation and preservation modes of Middle Pleistocene braided-river deposits were studied here based on in situ cosmogenic radionuclide concentrations. To account for potential discontinuous aggradation, a numerical model was developed to simulate the accumulation of cosmogenic radionuclides, 10Be and 26Al, in a sedimentary sequence and to account for deposition and erosion phases and post-depositional exposure. The method was applied to the Zutendaal gravels outcropping in NE Belgium, and 17 sediment samples were taken over a depth of 7 m and processed for determination of 10Be and 26Al concentrations. The model parameters were optimised using reduced χ2 minimisation. The Zutendaal gravels were deposited during (at least) three superimposed aggradational phases that were interrupted by stability or erosion lasting  40 kyr. This illustrates how long periods of non-deposition alternate with rapid and short depositional events. The key chronostratigraphical outcomes of this study, based on the optimal model outcomes, are as follows: (1) the fluvial deposits found at 7 m depth in As are dated at 654-62+218 ka, thereby corresponding to MIS 16; (2) the uppermost sequence of these deposits is dated at 562-45+211 ka, corresponding to MIS 14 and thereby possibly deposited by another stream apart from the Meuse River, which is assumed to have shifted its course east by that time; and (3) the total formation time of the upper 7 m of the Zutendaal gravels in As is about 90 kyr, so the deposition of the entire gravel sheet could correspond to more than one glacial period. The total formation time of the Zutendaal gravels constitutes nearly 20 % of the deposition age and shows the importance of considering the aggradation mode of braided rivers when applying CRN techniques to > 3 m deep sedimentary sequences.

Code availability

Python codes of the model, the optimisation scheme and the Gaussian kernel plot used in this study are open and available at (last access: 8 December 2022; Vandermaelen et al., 2022b).

Data availability

All data and code used in this paper are archived and publicly available in the UCLouvain DataVerse at (Vandermaelen et al., 2022c).

Author contributions

NV, VV and KB conceived the idea of applying CRN depth profiles to study the depositional history of the Campine Plateau. NV, VV and FC led the development of the numerical model. NV, KB, MC and VV participated in data acquisition and processing. KB and VV supervised the work and acquired the necessary funding. NV wrote the original draft, and all authors contributed to writing and revising the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors would like to thank Marco Bravin for assistance with manipulating and processing samples in the ELIc laboratory and Elise Dehaen and Rose Paque for the support during sample collection and UAV flights over the area. Furthermore, we thank Gilles Rixhon; an anonymous reviewer; and the handling associate editor, Hella Wittmann-Oelze, for their constructive comments on earlier versions of this paper.

Financial support

This research has been supported by a teaching assistantship grant from the Faculty of Sciences, UCLouvain, to Nathan Vandermaelen and a postdoctoral grant from the Fonds de la Recherche Scientifique (FRS-FNRS, Belgium) to François Clapuyt. This study was undertaken in the framework a research agreement between UCLouvain and SCK CEN (no. CO-19-17-4420-00; Belgian Nuclear Research Centre).

Review statement

This paper was edited by Hella Wittmann-Oelze and reviewed by Gilles Rixhon and one anonymous referee.


Akçar, N., Ivy-Ochs, S., Alfimov, V., Schlunegger, F., Claude, A., Reber, R., Christl, M., Vockenhuber, C., Dehnert, A., Rahn, M., and Schlüchter, C.: Isochron-burial dating of glaciofluvial deposits: First results from the Swiss Alps, Earth Surf. Process. Landf., 42, 2414–2425,, 2017. 

Balco, G. and Rovey, C. W.: An isochron method for cosmogenic-nuclide dating of buried soils and sediments, Am. J. Sci., 308, 1083–1114,, 2008. 

Balco, G., Stone, J. O. H., and Mason, J. A.: Numerical ages for Plio-Pleistocene glacial sediment sequences by 26Al /10Be dating of quartz in buried paleosols, Earth Planet. Sci. Lett., 232, 179–191,, 2005. 

Bats, H., Paulissen, E., and Jacobs, P.: De grindgroeve Hermans te As. Een beschermd landschap, Monumenten en Landschappen, 14, 56–63, 1995. 

Beerten, K., De Craen, M., and Wouters, L.: Patterns and estimates of post-Rupelian burial and erosion in the Campine area, north-eastern Belgium, Phys. Chem. Earth, 64, 12–20,, 2013. 

Beerten, K., Heyvaert, V. M. A., Vandenberghe, D., Van Nieuland, J., and Bogemans, F.: Revising the Gent Formation: a new lithostratigraphy for Quaternary wind-dominated sand deposits in Belgium, Geol. Belg., 20, 95–102,, 2017. 

Beerten, K., Dreesen, R., Janssen, J., and Van Uyten, D.: The Campine Plateau, in: Landscapes and Landforms of Belgium and Luxembourg, edited by: Demoulin, A., Springer, Berlin, Germany, 193–214,, 2018. 

Braucher, R., del Castillo, P., Siame, L., Hidy, A. J., and Bourlés, D. L.: Determination of both exposure time and erosion rate from an in situ-produced 10Be depth profile: A mathematical proof of uniqueness. Model sensitivity and applications to natural cases, Quat. Geochronol., 4, 56–67,, 2009. 

Braucher, R., Merchel, S., Borgomano, J., and Bourlès, D. L.: Production of cosmogenic radionuclides at great depth: A multi element approach, Earth Planet. Sci. Lett., 309, 1–9,, 2011. 

Bristow, C. S. and Best, J. L.: Braided rivers: perspectives and problems, in: Braided Rivers, Geological Society Special Publication No. 75, edited by: Best, J. L. and Bristow, C. S., Cambridge University Press, London, UK, 1-H,, 1993. 

Busschers, F. S., Kasse, C., van Balen, R. T., Vandenberghe, J., Cohen, K. M., Weerts, H. J. T., Wallinga, J., Johns, C., Cleveringa, P., and Bunnik, F. P. M.: Late Pleistocene evolution of the Rhine-Meuse system in the southern North Sea basin: imprints of climate change, sea-level oscillation and glacio-isostacy, Quat. Sci. Rev., 26, 3216–3248,, 2007. 

Chmeleff, J., von Blanckenburg, F., Kossert, K., and Jakob, D.: Determination of the 10Be half-life by multicollector ICP-MS and liquid scintillation counting, Nucl. Instrum. Methods Phys. Res. B: Beam Interact. Mater. At., 268, 192–199,, 2010. 

Christl, M., Vockenhuber, C., Kubik, P. W., Wacker, L., Lachner, J., Alfimov, V., and Synal, H. A.: The ETH Zurich AMS facilities: Performance parameters and reference materials, Nucl. Instrum. Methods Phys. Res. B: Beam Interact. Mater. At., 294, 29–38,, 2013. 

Claude, A., Akçar, N., Ivy-Ochs, S., Schlunegger, F., Kubik, P., Dehnert, A., Kuhlemann, J., Rahn, M., and Schlüchter, C.: Timing of early Quaternary accumulation in the Swiss Alpine Foreland, Geomorphology, 276, 71–85,, 2017. 

Cohen, K. M. and Gibbard, P.: Global chronostratigraphical correlation table for the last 2.7 million years, Subcommission on Quaternary Stratigraphy (International Commission on Stratigraphy), Cambridge, United Kingdom,, 2019. 

Covault, J. A., Craddock, W. H., Romans, B. W., Fildani, A., and Gosai, M.: Spatial and temporal variations in landscape evolution: Historic and longer-term sediment flux through global catchments, J. Geol., 121, 35–56,, 2013. 

De Brue, H., Poesen, J., and Notebaert, B.: What was the transport mode of large boulders in the Campine Plateau and the lower Meuse valley during the mid-Pleistocene?, Geomorphology, 228, 568–578,, 2015. 

Deckers, J., Vernes, R. W., Dabekaussen, W., Den Dulk, M., Doornenbal, J. C., Dusar, M., Hummelman, H. J., Matthijs, J., Menkovic, A., Reindersma, R. N., Walstra, J., Westerhoff, W. E., and Witmans, N.: Geologisch en hydrogeologisch 3D model van het Cenozoïcum van de Roerdalslenk in Zuidoost-Nederland en Vlaanderen (H3O-Roerdalslenk), Studie uitgevoerd door VITO, TNO-Geologische Dienst Nederland en de Belgische Geologische Dienst in opdracht van de Afdeling Land en Bodembescherming, Ondergrond, Natuurlijke Rijkdommen van de Vlaamse overheid, de Afdeling Operationeel Waterbeheer van de Vlaamse Milieumaatschappij, de Nederlandse Provincie Limburg en de Nederlandse Provincie Noord-Brabant, VITO en TNO-Geologische Dienst Nederland, VITO-rapport 2014/ETE/R/1, 205 pp., (last access: 24 February 2022), 2014. 

Dehaen, E.: Unraveling the characteristics of the Early and Middle Pleistocene Meuse River: study of the Zutendaal gravels on the Campine Plateau, MSc thesis, Faculty of Sciences, UCLouvain, Belgium, 63 pp., (last access: 28 February 2022), 2021. 

Dehnert, A., Kracht, O., Preusser, F., Akçar, N., Kemna, H. A., Kubik, P. W., and Schlüchter, C.: Cosmogenic isotope burial dating of fluvial sediments from the Lower Rhine Embayment, Germany, Quat. Geochronol., 6, 313–325,, 2011. 

Derese, C., Vandenberghe, D., Paulissen, E., and Van den haute, P.: Revisiting a type locality for Late Glacial aeolian sand deposition in NW Europe: Optical dating of the dune complex at Opgrimbie (NE Belgium), Geomorphology, 109, 27–35,, 2009. 

Dunai, T. J. (Ed.): Cosmogenic Nuclides, Cambridge University Press, New York, USA,, 2010. 

Erlanger, E. D., Granger, D. E., and Gibbon, R. J.: Rock uplift rates in South Africa from isochron burial dating of fluvial and marine terraces, Geology, 40, 1019–1022,, 2012. 

GDI-Vlaanderen: Digitaal Hoogtemodel Vlaanderen II, Version 2014.01,, last access: 24 February 2022. 

Granger, D. E. and Muzikar, P. F.: Dating sediment burial with in situ-produced cosmogenic nuclides: theory, techniques, and limitations, Earth Planet. Sci. Lett., 188, 269–281,, 2001. 

Gullentops, F., Janssen, J., and Paulissen, E.: Saalian nivation activity in the Bosbeek valley, NE Belgium, Geologie en Mijnbouw, 72, 125–130, 1993. 

Gullentops, F., Bogemans, F., de Moor, G., Paulissen, E., and Pissart, A.: Quaternary lithostratigraphic units (Belgium), Geol. Belg., 4, 153–164,, 2001. 

Hancock, G. S., Anderson, R. S., Chadwick, O. A., and Finkel, R. C.: Dating fluvial terraces with 10Be and 26Al profiles: application to the Wind River, Wyoming, Geomorphology, 27, 41–60,, 1999. 

Hidy, A. J., Gosse, J. C., Pederson, J. L., Mattern, J. P., and Finkel, R. C.: A geologically constrained Monte Carlo approach to modeling exposure ages from profiles of cosmogenic nuclides: An example from Lees Ferry, Arizona, Geom. Geophys., 11, Q0AA10,, 2010. 

Hidy, A. J., Gosse, J. C., Sanborn, P., and Froese, D. G.: Age-erosion constraints on an Early Pleistocene paleosol in Yukon, Canada, with profiles of 10Be and 26Al: Evidence for a significant loess cover effect on cosmogenic nuclide production rates, Catena, 165, 260–271,, 2018. 

Knudsen, M. F., Egholm, D. L., and Jansen, J. D.: Time-integrating cosmogenic nuclide inventories under the influence of variable erosion, exposure, and sediment mixing, Quat. Geochronol. 51, 110–119,, 2019. 

Korschinek, G., Bergmaier, A., Faestermann, T., Gerstmann, U. C., Knie, K., Rugel, G., Wallner, A., Dillmann, I., Dollinger, G., Lierse Von Gostomski, C., Kossert, K., Maiti, M., Poutivtsev, M., and Remmert, A.: A new value for the half-life of 10Be by Heavy-Ion Elastic Recoil Detection and liquid scintillation counting, Nucl. Instrum. Methods Phys. Res. B: Beam Interact. Mater. At., 268, 187–191,, 2010. 

Lachner, J., Christl, M., Müller, A. M., Suter, M., and Synal, H. A.: 10Be and 26Al low-energy AMS using He-stripping and background suppression via an absorber, Nucl. Instrum. Methods Phys. Res. B: Beam Interact. Mater. At., 331, 209–214,, 2014. 

Laloy, E., Beerten, K., Vanacker, V., Christl, M., Rogiers, B., and Wouters, L.: Bayesian inversion of a CRN depth profile to infer Quaternary erosion of the northwestern Campine Plateau (NE Belgium), Earth Surf. Dynam., 5, 331–345,, 2017. 

Lauer, T., Frechen, M., Hoselmann, C., and Tsukamoto, S.: Fluvial aggradation phases in the Upper Rhine Graben-new insights by quartz OSL dating, Proc. Geol. Assoc., 121, 154–161,, 2010. 

Lauer, T., Weiss, M., Bernhardt, W., Heinrich, S., Rappsilber, I., Stahlschmidt, M. C., von Suchodoletz, H., and Wansa, S.: The Middle Pleistocene fluvial sequence at Uichteritz, central Germany: Chronological framework, paleoenvironmental history and early human presence during MIS 11, Geomorphology, 354, 107016,, 2020. 

Le Dortz, K., Meyer, B., Sébrier, M., Braucher, R., Nazari, H., Benedetti, L., Fattahi, M., Bourlès, D., Foroutan, M., Siame, L., Rashidi, A., and Bateman, M. D.: Dating inset terraces and offset fans along the Dehshir Fault (Iran) combining cosmogenic and OSL methods, Geophys. J. Int., 185, 1147–1174,, 2011. 

Margreth, A., Gosse, J. C., and Dyke, A. S.: Quantification of subaerial and episodic subglacial erosion rates on high latitude upland plateaus: Cumberland Peninsula, Baffin Island, Arctic Canada, Quat. Sci. Rev. 133, 108–129,, 2016. 

Martin, L. C. P., Blard, P. H., Balco, G., Lavé, J., Delunel, R., Lifton, N., and Laurent, V.: The CREp program and the ICE-D production rate calibration database: A fully parameterizable and updated online tool to compute cosmic-ray exposure ages, Quat. Geochronol., 38, 25–49,, 2017. 

Miall, A. D. (Ed.): The Geology of fluvial deposits, Springer, Berlin, Germany,, 1996. 

Mol, J., Vandenberghe, J., and Kasse, C.: River response to variations of periglacial climate in mid-latitude Europe, Geomorphology, 33, 131–148,, 2000. 

Nichols, K. K., Bierman, P. R., Hooke, R. L., Clapp, E. M., and Caffee, M.: Quantifying sediment transport on desert piedmonts using 10Be and 26Al, Geomorphology, 45, 105–125, 2002. 

Nichols, K. K., Bierman, P. R., Eppes, M. C., Caffee, M., Finkel, R., and Larsen, J.: Late Quaternary history of the Chemehuevi Mountain Piedmont, Mojave Desert, deciphered using 10Be and 26Al, Am. J. Sci., 305, 345–368,, 2005. 

Nishiizumi, K.: Cosmic ray production rates of 10Be and 26Al in quartz from glacially polished rocks, J. Geophys. Res.-Sol. Ea., 94, 17907–17915,, 1989. 

Nishiizumi, K.: Preparation of 26Al AMS standards, J Nucl. Instrum. Methods Phys. Res. B: Beam Interact. Mater. At., 223–224, 388–392,, 2004. 

Pannekoek, A. J.: Einigen Notizen über die Terrassen in Mittel- und Nord-Limburg, Natuurhistorisch Maandblad, 13, 89–92, 1924. 

Paulissen, E.: De Morfologie en de Kwartairstratigrafie van de Maasvallei in Belgisch Limburg, Verhandelingen van de koninklijke Vlaamse academie voorwetenschappen, letteren en schone kunsten van België, Klasse der Wetenschappen, 127, 1–266, 1973. 

Paulissen, E.: Les nappes alluviales et les failles Quaternaires du Plateau de Campine, in: Guides Géologiques Régionaux – Belgique, edited by: Robaszynski, F. and Dupuis, C., Masson, Paris, France, 167–170, ISBN10 2225753946, ISBN13 9782225753947, 1983. 

Portenga, E. W. and Bierman, P. R.: Understanding earth's eroding surface with 10Be, GSA Today, 21, 4–10,, 2011. 

Rixhon, G., Braucher, R., Bourlès, D., Siame, L., Bovy, B., and Demoulin, A.: Quaternary river incision in NE Ardennes (Belgium) – Insights from 10Be /26Al dating of river terraces, Quat. Geochronol., 6, 273–284,, 2011. 

Rixhon, G., Bourlès, D. L., Braucher, R., Siame, L., Cordy, J. M., and Demoulin, A.: 10Be dating of the Main Terrace level in the Amblève valley (Ardennes, Belgium): New age constraint on the archaeological and palaeontological filling of the Belle-Roche palaeokarst, Boreas, 43, 528–542,, 2014. 

Rizza, M., Abdrakhmatov, K., Walker, R., Braucher, R., Guillou, V., Carr, A. S., Campbell, G., McKenzie, D., Jackson, J., Aumaître, G., Bourlès, D. L., and Keddadouche, K.: Rate of slip from multiple Quaternary dating methods and paleoseismic investigations along the Talas-Fergana Fault: tectonic implications for the Tien Shan Range, Tectonics, 38, 2477–2505,, 2019. 

Rodés, A., Pallàs, R., Braucher, R., Moreno, X., Masana, E., and Bourlès, D.: Effect of density uncertainties in cosmogenic 10Be depth-profiles: Dating a cemented Pleistocene alluvial fan (Carboneras Fault, SE Iberia), Quat. Geochronol., 6, 186–194,, 2011. 

Schaller, M., von Blanckenburg, F., Hovius, N., and Kubik, P. W.: Large-scale erosion rates from in situ-produced cosmogenic nuclides in European river sediments, Earth Planet. Sci. Lett., 188, 441–458,, 2001. 

Schaller, M., Ehlers, T. A., Blum, J. D., and Kallenberg, M. A.: Quantifying glacial moraine age, denudation, and soil mixing with cosmogenic nuclide depth profiles, J. Geophys. Res. 114, F01012,, 2009. 

Schoonejans, J., Vanacker, V., Opfergelt, S., Granet, M., and Chabaux, F.: Coupling uranium series and 10Be cosmogenic radionuclides to evaluate steady-state soil thickness in the Betic Cordillera, Chem. Geol., 446, 99–109,, 2016. 

Siame, L., Bellier, O., Braucher, R., Sébrier, M., Cushing, M., Bourlès, D., Hamelin, B., Baroux, E., de Voogd, B., Raisbeck, G., and Yiou, F.: Local erosion rates versus active tectonics: Cosmic ray exposure modelling in Provence (south-east France), Earth Planet. Sci. Lett., 220, 345–364,, 2004. 

Sougnez, N. and Vanacker, V.: The topographic signature of Quaternary tectonic uplift in the Ardennes massif (Western Europe), Hydrol. Earth Syst. Sci., 15, 1095–1107,, 2011. 

Stone, J. O.: Air pressure and cosmogenic isotope production, J. Geophys. Res.-Sol. Ea., 105, 23753–23759,, 2000. 

Taylor, J. R. (Ed.): An introduction to error analysis, University science books, Sausalito, California, USA, ISBN10 9780935702750, ISBN13 978-0935702750, 1997. 

U.S. Geological Survey's EROS Data Center: Global 30 Arc-Second Elevation (GTOPO30),, 1996. 

Vanacker, V., von Blanckenburg, F., Hewawasam, T., and Kubik, P. W.: Constraining landscape development of the Sri Lankan escarpment with cosmogenic nuclides in river sediment, Earth Planet. Sci. Lett., 253, 402–414,, 2007. 

Vanacker, V., von Blanckenburg, F., Govers, G., Molina, A., Campforts, B., and Kubik, P. W.: Transient river response, captured by channel steepness and its concavity, Geomorphology, 228, 234–243,, 2015. 

Van Balen, R. T., Houtgast, R. F., van der Wateren, F. M., Vandenberghe, J., and Bogaart, P. W.: Sediment budget and tectonic evolution of the Meuse catchment in the Ardennes and the Roer Valley Rift System, Glob. Planet Change​​​​​​​, 27, 113–129,, 2000. 

van den Berg, M.: Fluvial sequences of the Meuse – a 10 Ma record of neotectonics and climate change at various time-scales, PhD thesis, Wageningen University, 181 pp., (last access: 6 March 2022), 1996. 

Vandenberghe, D., Vanneste, K., Verbeeck, K., Paulissen, E., Buylaert, J.-P., De Corte, F., and Van den haute, P.: Late Weichselian and Holocene earthquake events along the Geleen fault in NE Belgium: OSL age constraints, Quat. Int., 199, 56–74,, 2009. 

Vandenberghe, J.: Timescale, Climate and River Development, Quat. Sci. Rev., 14, 631–639,​​​​​​​, 1995. 

Vandenberghe, J.: A typology of Pleistocene cold-based rivers, Quat. Int., 79, 111–121,​​​​​​​, 2001. 

Vandermaelen, N., Vanacker, V., Clapuyt, F., Christl, M., and Beerten, K.: Reconstructing the depositional history of Pleistocene fluvial deposits based on grain size, elemental geochemistry and in situ 10Be data, Geomorphology, 402, 108127,, 2022a. 

Vandermaelen, N., Clapuyt, F., and Vanacker, V.: Repository hosting resources for numerical modelling of complex CRN accumulation history, Github [code], available at:, last access: 8 December 2022b. 

Vandermaelen, N., Beerten, K., Clapuyt, F., Christl, M., and Vanacker, V.: CRN Datasets and numerical modeling tools for geochronological study at geosite Hermans on the Campine Plateau, UCLouvain DataVerse [data set],, last access: 8 December 2022c. 

Vanneste, K., Verbeeck, K., Camelbeeck, T., Paulissen, E., Meghraoui, M., Renardy, F., Jongmans, D., and Frechen, M.: Surface-rupturing history of the Bree fault scarp,Roer Valley graben: Evidence for six events since the late Pleistocene, J. Seismol., 5, 329–359, 2001. 

von Blanckenburg, F., Belshaw, N. S., and O'Nions, R. K.: Separation of 9Be and cosmogenic 10Be from environmental materials and SIMS isotope dilution analysis, Chem. Geol., 129, 93–99,, 1996. 

Westerhoff, W. E., Kemna, H. A., and Boenigk, W.: The confluence area of Rhine, Meuse, and Belgian rivers: Late Pliocene and Early Pleistocene fluvial history of the northern Lower Rhine Embayment, Neth. J. Geosci., 87, 107–125,, 2008.  

Xu, L., Ran, Y., Liu, H., and Li, A.: 10Be-derived sub-Milankovitch chronology of Late Pleistocene alluvial terraces along the piedmont of SW Tian Shan, Geomorphology, 328, 173–182,, 2019. 

Short summary
We constrained deposition phases of fluvial sediments (NE Belgium) over the last 1 Myr with analysis and modelling of rare isotopes accumulation within sediments, occurring as a function of time and inverse function of depth. They allowed the determination of three superposed deposition phases and intercalated non-deposition periods of ~ 40 kyr each. These phases correspond to 20 % of the sediment age, which highlights the importance of considering deposition phase when dating fluvial sediments.