A new model for radiocarbon dating of marine shells from the Netherlands

Abstract

A new model for the interpretation of radiocarbon (¹⁴C) dates of Holocene marine shells is presented. For the Netherlands, the size of reservoir effect is difficult to assess, as these shells often lived in an environment of mixed marine- and river waters. Both stable isotopes ¹³C and ¹⁸O of the shell carbonate give insight in the environmental conditions the shells lived in. River water occurs in two main categories, distinguished by ¹⁸O: the Rhine which is dominant, and other rivers. This leads to two estuary mixing lines between the North Sea and rivers. The stable isotopes of the shell carbonate are also indicative for additional processes, such as uptake of secondary carbonate from the soil by shells, and exchange of C isotopes between atmosphere and water. Extensions of the main model deal with special cases such as pools of stagnant water and lakes. The model leads to an assessment of the recent ¹⁴C activities of the system the shells lived in, called ¹⁴a^SYS. The measured ¹⁴C activities relative to these ¹⁴a^SYS values determine the ¹⁴C age of the shells and include the reservoir effect. This way we circumvent normalizing to δ¹³C = –25‰, i.e. the terrestrial timescale and the subsequent correction for reservoir effects. The model is applied to a large legacy dataset of marine shells from the Netherlands, obtained during the last 7 decades. It contains 1116 ¹⁴C dated shells; for the majority of these, the 3 isotopes ¹³C, ¹⁴C and ¹⁸O are measured.

Introduction

Radiocarbon (¹⁴C) dates are by convention normalized for fractionation using the stable carbon isotope to δ¹³C = –25‰. Together with other normalizations (concerning half-life and standard activity), this defines the timescale in BP (Mook and van der Plicht 1999; Stuiver and Polach 1977). This timescale is valid for terrestrial samples. Aquatic shells are subject to reservoir effects. Generally, aquatic reservoirs contain less ¹⁴C than terrestrial reservoirs and the atmosphere, causing an offset in ¹⁴C dates. Therefore, when aquatic dates are normalized for fractionation and expressed in BP, subtraction of reservoir effects is needed for a correct age estimation (e.g. Cook et al. 2001; Lanting and van der Plicht 1998; Olsson 1983; Philippsen 2013).

For marine surface samples, a general value of 400 years applies for the size of the reservoir effect for the Southern North Sea during the Holocene. Our dataset concerns this locality and timerange. For other regions and time periods, the magnitude of the reservoir effect can be different (Alves et al. 2018). We note that for the present calibration curve Marine20 (which is model based) it is more like 550 years (Heaton et al. 2020, 2023).

For absolute shell dates (BC, AD or calBP) this 400 years needs to be subtracted (from the measured date in BP) before calibrating using the terrestrial calibration curve IntCal20 (Reimer et al. 2020). Alternatively, one can calibrate the BP date using the marine calibration curve Marine20 (Heaton et al. 2020), which incorporates the reservoir effect. Marine20 takes into account temporal and local offsets from the global average value. The marine reservoir effect is caused by upwelling of deep ocean waters which are old and thus depleted in ¹⁴C. The surface ocean exchanges CO₂ (including ¹⁴CO₂) with the atmosphere but is not recent (in terms of ¹⁴C) because of mixture with this upwelling deep and old ocean water.

The general “400 years rule” only works for pure marine shells. For other aquatic reservoirs like rivers different values apply. This reservoir effect is mainly caused by dissolved fossil (¹⁴C-free) carbonate. For the Netherlands, the dominant river is the Rhine which shows reservoir effects of around 1300 years (Lanting and van der Plicht 1998; Mook 2006). Other water bodies may show different values. This is the case for estuaries and lagoons, which consist of a mixture of sea- and river waters.

There is no general numerical reservoir effect value that can be used for shells from aquatic waters other than 100% marine, such as rivers, lakes, lagoons, estuaries and other mixed water bodies. Only in exceptional cases one can compare ¹⁴C dates of shells directly with those from terrestrial samples. Examples are shells from a well defined archaeological context or geological stratum. There are even cases where there is no (or small) reservoir effect, such as stagnant waters which can reach CO₂ equilibrium with the atmosphere (e.g. van der Plicht et al. 2020). Hence, transformation of ¹⁴C dates of shells from not-100% marine environments into absolute dates is problematic. This is the case for several hunderds of shell dates measured over the years (since 1952) in Groningen. The dataset consists of marine shells, as well as many shells from rivers, lakes, stagnant waters, estuaries and a variety of mixed water bodies.

We developed a model for obtaining absolute dates for such shells. It uses both stable isotopes ¹³C and ¹⁸O of the shell carbonate as tracers to reconstruct the isotopic environment where the shell species lived in, taking into account mixing of waters, fractionation and isotopic exchange processes. This way we can circumvent the need for reservoir effect correction, in fact by not normalizing to the terrestrial timescale. This more or less follows the practice used in isotope hydrology (Clark and Fritz 1997; Mook 2006). The modeled ¹⁴C activities result in dates that can be made absolute by using the calibration curve IntCal20 (Reimer et al. 2020). These modeled dates are equivalent to but should not be expressed in the unit BP because they are not normalized for fractionation; they are formally not conventional. We will use here “¹⁴C years” or “modeled years.”

The present work is enabled by our large dataset of shells. It is unique by virtue of its size, and a large amount of ¹⁸O measurements. The latter isotope was not commonly measured in the first decades of radiocarbon dating. Of importance is also that the shells originate from a well studied delta including much archaeological and geophysical knowledge for a relatively small region. The dataset provides proper testing of our model.

The large density of measurements even allows to estimate the δ¹⁸O value for the shells for which only the δ¹³C value has been measured. This enables the application of our model to derive ¹⁴C ages (albeit with increased uncertainty). The full dataset is available as supporting online material. It includes all shells analyzed in Groningen during the years 1952–2023. Detailed applications concerning past landscape developments in the Netherlands will be published in the more specialized or local geophysical/archaeological literature.

Dataset and methods

The shell carbonate dataset measured in Groningen contains 1116 shells of the Netherlands which are dated by ¹⁴C. For about 70% of these, the 3 isotopes ¹³C, ¹⁴C and ¹⁸O are measured.

There is a significant amount of shells for which only the two carbon isotopes are measured.

In addition, 20 shell dates measured in the 1950s do not have stable isotope data at all, also not ¹³C. These data are valuable nevertheless, the lacking stable isotopes can often be estimated (justified by the mere size of the total dataset) and subsequently be analyzed by our model.

In this article we focus on shells with the 3 isotopes available from a Holocene geophysical context, which amounts to about 700 shells. The remaining data concern shells from an archaeological context, recent shells and a few shells which appeared to be of Pleistocene age.

Over the years, a selection of measurements found their way in publications and theses. These were usually interpreted in the users community using a reservoir effect of 400 years. For completeness, also included is a dataset of 264 shells for which only the stable isotopes ¹³C and ¹⁸O are measured; these were mostly analyzed before the introduction of AMS.

A modern map of the Netherlands is shown in Figure 1. The datatables are organized according to regions. These regions correspond to provinces, except for the Wadden Islands, which belong to three provinces and are taken together here. Figure 1 also shows the main water bodies: the North Sea and the rivers Rhine, Meuse and Scheldt. The Rhine is the dominant river; the flow rate is about 2000 m³/s. For the Meuse, that is 230 m³/s; for the Scheldt, 120 m³/s. These numbers indicate the annual averages. The IJssel river is a separate branch of the Rhine. The Lake IJssel is a freshwater inner lake; it used to be part of the North Sea called Zuiderzee. Since 1932 it is disconnected from the sea by a 32 km long closure dike (called Afsluitdijk). The former lake area is partly used for land reclamation and is today the province of Flevoland.

Figure 1. Modern map of the Netherlands, showing the main water bodies: the North Sea, Waddensea and Lake IJssel, and the rivers Scheldt, Meuse, Rhine and IJssel. The coastal provinces are indicated, as well as the closure dike (Afsluitdijk) between North Sea/Waddensea and Lake IJssel, the former Zuiderzee. “A” is the location of Amsterdam.

For illustrative purposes, a map showing a reconstruction of the country dating 3850 BC is included (Figure 2). The map is taken from Vos (2015). Maps from other times during the Holocene are also available, these can be obtained from the Rijksdienst voor het Cultureel Erfgoed (2024). Together they provide an excellent illustration of the temporal changes in the prehistoric landscape during the Holocene, in particular the development of coastline, lakes, estuaries and course of rivers.

Figure 2. Map of the Netherlands dating ca. 3850 BC. The modern provinces are shown by the black lines. The map is taken from Vos (2015), which reference also contain maps from other times during the Holocene. The rivers Rhine, Meuse, Scheldt, Vecht and Eems are indicated.

The ¹⁴C shell dates were obtained by 4 different instrumental methodologies: the conventional method (laboratory code GrN, used until 2011; Mook and Streurman 1983), AMS using the Tandetron (laboratory code GrA, used 1994-2017; Aerts et al. 2001; van der Plicht et al. 2000) and the Micadas (laboratory code GrM, used since 2017; Dee et al. 2020). Before Groningen had an operational AMS facility, shells were pretreated and graphitized in Groningen and subsequently AMS-dated in Uppsala (laboratory code Ua; van Klinken 1991). The stable isotope ratios were all measured by Isotope Ratio Mass Spectrometry (IRMS) in Groningen.

The stable isotope ratios are expressed in a delta value, which for ¹³C is a measure of the difference of the ¹³C/¹²C isotope ratio of a sample and that of a reference material:

([1]) $$\rm {\rm ^{13}\delta = {\delta^{13}}C = {\left( {^{13}C{/^{12}}C} \right)_{sample}}/ \ {\left( {^{13}C{/^{12}}C} \right)_{reference}}- 1 ( \times 1000 \unicode{x2030}))}$$

Similar, for ¹⁸O this is:

([2]) $$\rm{ ^{18}\delta = {\delta^{18}}O = {\left( {^{18}O{/^{16}}O} \right)_{sample}}/ \ {\left( {^{18}O{/^{16}}O} \right)_{reference}}- 1 ( \times 1000 \unicode{x2030})}$$

For both ¹³C and ¹⁸O in carbonates, the reference is an international standard limestone known as Pee Dee Belemnite, or PDB which has an accurately measured absolute ¹³C concentration. For ¹⁸O in water, the reference is SMOW (standard mean ocean water). The reference materials have accurately measured absolute isotope concentrations. For more details, see Mook (2006) and references therein. The short notations ¹³δ and ¹⁸δ will be used in the more complex formulas in order to make them better readable. The measurement precision of the stable isotope ratios is about 0.1‰.

The ¹⁴C measurements are reported in activity ratios ¹⁴a, as defined by Mook and van der Plicht (1999). It follows the convention by using the half-life value of 5568 years (the Libby value) and standard activity (95% of the activity of the specific batch of Oxalic Acid nr.1 in 1950).

The convention also requires a correction for isotopic fractionation to δ¹³C = −25‰, here called normalization using

([3]) $$\rm{ ^{14}{a_N}{\ =\ ^{14}}a [1-2({\delta^{13}}C + 0.025)]}$$

which is valid as a good approximation (Mook 2006). The δ¹³C reference value of −25‰ is chosen as typical for charcoal, wood and peat of C3 plants (Mook and Streurman 1983). The normalized activity ratio ¹⁴a_N is equivalent to F¹⁴C; see e.g. van der Plicht and Hogg (2006) for the nomenclature. Both lead to ¹⁴C dates in BP, the terrestrial timescale.

This normalization leads to an offset with dates of aquatic origin, known as the reservoir effect. For the model introduced in this manuscript, we use not-normalized ¹⁴C activity ratios, ¹⁴a (and not ¹⁴a_N) in order to avoid the need for reservoir corrections.

For conventional samples, the ¹⁴C activities are directly measured in CO₂ by proportional gas counters, leading to ¹⁴a values (using the conventional half-life and oxalic acid standard). For the same gas, also the δ¹³C value is measured by IRMS (Mook and Streurman 1983).

For AMS samples however, ambiguities do exist. The CO₂ gas is transferred into graphite; the AMS measures the ¹³C and ¹⁴C concentration in the graphite (Dee et al. 2020). The AMS dates are subsequently normalized using the δ¹³C value as measured by the AMS. This takes into account the fractionation caused by the graphite production and in the accelerator system, which can be significant.

The δ¹³C value for all samples is also measured by IRMS. Hence there are two measurements available for the δ¹³C values. The AMS δ¹³C value is less precise than the one measured by IRMS, typically by about one order of magnitude. The IRMS value is reported, the one measured by AMS is not reported. The IRMS δ¹³C value represents the proper ¹³C value of the shell; the AMS ¹⁴a_N value or BP date is calculated using the AMS δ¹³C value. For normalized dates this is not an issue. However, our model is based on ¹⁴a (i.e. not ¹⁴a_N). The problem is solved by de-normalizing the AMS ¹⁴a_N value using the IRMS-measured proper δ¹³C value:

([4]) $$\rm{^{14}a {\ =\ ^{14}}{a_N}{,_{AMS}}/\ [1-2({\delta^{13}}{C_{IRMS}} + 0.025)]}$$

The effect is not large, usually well within the measurement uncertainty.

The modeled ¹⁴C dates, based on not-normalized ¹⁴C activities are calibrated using IntCal20 (Reimer et al. 2020) and the OxCal code (version 4.4, https://c14.arch.ox.ac.uk/oxcal/OxCal.html). In our tables the calibrated dates are reported in calBP, i.e. calender years relative to 1950 AD. Reported are the 95.4% probability date ranges, rounded to 5 years and taking separated date ranges together.

Isotopes in water and carbonate

Our model is based on the isotopes ¹³C, ¹⁴C and ¹⁸O in natural waters and carbonate. The general isotope values for the main reservoirs relevant for the Netherlands are shown in Table 1. The table distinguishes the three main water bodies marine (the North Sea) and rivers (Rhine and others).

Table 1. General contents of the stable isotopes ¹³C and ¹⁸O for water and shell carbonate in the Netherlands. For water, the values of Chlorinity [Cl] and modern ¹⁴C activity (named ¹⁴a^SYS, see text) are also shown

The chlorinity of the water is shown as [Cl]. This [Cl] value is defined as the concentration of chloride (plus bromide and iodide) given in g/kg or ‰ (Mook 1971, 2006). For the sea, the Cl concentration is 560 mmol/L, which corresponds to 19.36 ppt (g/kg) or 19.36‰. Freshwater does not contain Cl, hence its concentration is 0 mmol/L or 0‰.

For North Sea water, measurements show that for dissolved inorganic carbon (DIC) the typical value for δ¹³C is 1.5‰ which is based on measurements by Mook (1968, 104). The δ¹⁸O value is 0 (the SMOW reference value). Also, for shell carbonate, the δ¹³C values for a large body of North Sea shells have been measured, resulting in a typical value of 1.5‰. For δ¹⁸O, the fractionation between water and shell carbonate is known to be 1‰ (Mook 1968, 2006).

In terms of stable isotopes, the rivers in the Netherlands can be distinguished in two main categories: the Rhine, and others. The δ¹³C values for bicarbonate and shell carbonate for all freshwater bodies is –10.5‰. For ¹⁸O however, there is a difference between various rivers. The Rhine is dominant and is mainly fed by snow from the Alp mountains in Switzerland. The ¹⁸O contents is therefore depleted by the so-called altitude effect, which is basically a temperature effect. About 800 km downstream and including contributions from tributary rivers along the way, in the Netherlands the typical δ¹⁸O value for Rhinewater is –10.0‰ (SMOW).

The other rivers (Meuse, Scheldt and other smaller streams) are rain-fed by precipitation at low altitudes. Their average δ¹⁸O value for water is –7.3‰ (SMOW) (Mook 1970). The 1‰ fractionation between water and shell carbonate leads for shells to δ¹⁸O = −9.0‰ (PDB) for the Rhine, and −6.3‰ (PDB) for the other (precipitation fed) rivers. This forms the basis of what we call “model R” and “model P”, respectively. The values δ¹⁸O = –10.0‰ and –7.3‰ for water are measured annual averages. They show small seasonal variations with an amplitude of about 0.5‰ (Mook 1970, 2006).

The numbers shown in Table 1 are widely used in hydrological applications. For a more detailed treatise we refer to the literature (Clark and Fritz 1997; Mook 2006; see also IAEA 2024 and references therein). For radiocarbon, the typical value for the recent ¹⁴C activity ¹⁴a^SYS for river water observed throughout the Netherlands is 85%. A selection of measurements leading to this number is given in Appendix A. This recent activity corresponds to a riverine reservoir effect of about 1300 years.

The shell dating model

This chapter describes in detail a new ¹⁴C dating model for shells. It is based on the three isotopes ¹³C, ¹⁴C and ¹⁸O of the shell carbonate.

The concentrations of the stable isotopes ¹³C and ¹⁸O are different for shells living in different environments. There are marine shells, riverine shells and shells living in estuaries, consisting of mixed river and sea water. Riverwater can be distinguished in two categories: Rhine water, and water from other rivers. The stable isotope contents in shells also depends on additions effects, such as isotopic exchange between water and the atmosphere, and uptake of carbonate from the soil.

There are other aquatic environments which are habitats for shells: stagnant waters and lakes, and various mixed water bodies. The latter can also show variations through time which obviously complicates the dating procedure.

The isotopic behavior of water and shells is investigated for all the above mentioned environments, including equations for the three natural isotopes involved. This then leads to our dating model for shells, which includes handling reservoir effects.

Following a general introduction, the various environments and isotopic effects therein are discussed below in sections a) through g).

The basic ¹⁴C age equations, using the conventional half-life T_1/2 = 5568 years are the following:

([5]) $$\rm{ ^{14}a\left( t \right) {\ =\ ^{14}}a\left( {t = 0} \right)\ exp(-\lambda t)\ with\ \lambda = ln2/{T_{1/2}} = 0.693/5568}$$

([6]) $${\rm age (t) = -8033\ ln\left[ {^{14}a\left( t \right){/^{14}}a\left( {t = 0} \right)} \right]}$$

The time t = 0 represents the time of death of the organism, the “start of the ¹⁴C clock” in dating terms. This corresponds to the “modern” or “recent” activity, i.e. the ¹⁴C content of the environment the organism was living in. This in turn also represents (in principle) the activity of the standard, ¹⁴a(t = 0) = ¹⁴a⁰ (Mook and van der Plicht 1999), refering to “year 0,” i.e. 1950, the reference year.

When normalized for isotopic fractionation, the terrestrial recent value (t = 0) will remain the same: ¹⁴a_N = 100%. However, the marine value becomes ¹⁴a_N = 105%. The difference of 25‰ in δ¹³C corresponds to 50‰ = 5% in ¹⁴a or (approximately) 400 ¹⁴C years (for not-too-old samples).

We developed a new method for ¹⁴C dating for interpretation of dates from aquatic samples, instigated by our large series of shells from the Netherlands, which is in fact a large delta region (see Figure 1). These shells live in environments consisting of water bodies with various compositions. The marine water is North Sea water, for which the reservoir effect is 400 years during the Holocene. The freshwaters can be rivers, lakes or lagunes leading to complex isotopic mixing situations of water bodies, and thus unknown reservoir effects.

We introduce here the symbol ¹⁴a^SYS (SYS for system or reservoir) which is more consistent for our purpose, directly indicating the original natural ¹⁴C activity ratio of the system (environment) the organism to be dated lived in. This is not necessarily the same as ¹⁴a⁰ as defined in radiometry, or the number defined in the ¹⁴C convention. This way, for a recent terrestrial sample the reference values are then ¹⁴a^SYS = 100% and δ¹³C = −25‰. For a recent fully marine sample, these values are ¹⁴a^SYS = 100% and δ¹³C = 0‰.

Thus, our model is used for interpretation of ¹⁴C dates for aquatic samples without normalization to δ¹³C = −25‰, i.e. the terrestrial regime which makes subtracting reservoir effects unnecessary. Instead, we will use tracers to follow isotopic processes: in particular δ¹⁸O values of the dated shells, basic water properties such as Chlorinity, and specific shell species behavior.

The starting point will be the ¹⁴a^SYS value, which would have been the ¹⁴a value measured at the time the organism to be dated lived in.

For a sample with measured ¹⁴C activity ¹⁴a, the ¹⁴C age is then calculated as

([7]) $${\rm age = -8033ln\left( {^{14}a{/^{14}}{a^{SYS}}} \right)}$$

using the conventional half-life T_1/2 = 5568 years and the oxalic acid reference. The normalization to δ¹³C = –25‰ (fractionation correction) is not applied here. Again, this age is therefore formally not expressed in BP.

The formula [7] can also be written as:

([8]) $${\rm age = -8033ln\left( {^{14}a{/^{14}}{a^{SYS}}} \right) = \left[ {-8033ln\left( {^{14}a} \right)} \right]-\left[ {-8033ln\left( {^{14}{a^{SYS}}} \right)} \right] = -8033ln\left( {^{14}a} \right)-R\;}$$

where R = –8033ln(¹⁴a^SYS), which can be interpreted as a measure of the reservoir effect.

As an illustrative numerical example, the recent activity for freshwater rivers is ¹⁴a^SYS = 85% giving R = 1305 years. For the North Sea, ¹⁴a^SYS = 100% giving R = 0; this obviously follows from using the not normalized activity. The value for normalized sea water is ¹⁴a_N = 105%, which results in –8033ln(¹⁴a_N) = –392 years.

These numbers correspond (well within uncertainties) with the general rounded values of 1300 and 400 years used for the freshwater and marine reservoir effects, respectively.

a) Estuary mixing

For estuaries, the stable isotope values of the shell carbonate (numbers shown in Table 1) lead to the following mixing lines. These mixing equations are possible because the concentration of bicarbonate [HCO₃ ^–] for North Sea and rivers are the same, 2.5 mmol/L.

There are some exceptional cases, but the effect can be neglected within uncertainties.

The mixing relations for the stable isotopes are

([9]) $${\rm {\delta^{13}}C = 1.64\ {\delta^{18}}O-0.14\ for\ precipitation\ fed\ rivers\ \left( {model \,\,P} \right)}$$

based on the endpoints (1.0, 1.5) and (–6.3, –10.5) for (δ¹⁸O, δ¹³C); and

([10]) $${\rm {\delta^{13}}C = 1.20\ {\delta^{18}}O + 0.30\ for\ the\ Rhine\ river\ \left( {model\,\, R} \right)}$$

based on the endpoints (1.0, 1.5) and (–9.0, –10.5) for (δ¹⁸O, δ¹³C).

The relation between the ¹⁴C and ¹³C is another linear mixing relation, which can be expressed as follows:

([11]) $${\rm {(^{14}}{a^{SYS}}{-^{14}}{a_m})/{(^{14}}{a_f}{-^{14}}{a_m}) = {(^{13}}\delta{-^{13}}{\delta_m})/{(^{13}}{\delta_f}{-^{13}}{\delta_m}) = {(^{18}}\delta{-^{18}}{\delta_m})/{(^{18}}{\delta_f}{-^{18}}{\delta_m})}$$

where “m” stands for marine-, and “f” for freshwater. Using the numbers of Table 1, the C isotope mixing formula for shell carbonate is (with ¹⁴a in % and ¹³δ in ‰) as follows:

([12]) $$\rm{ ^{14}{a^{SYS}} = {1.25^{13}}\delta + 98.125\;\;\;\;\;\;}$$

An analogue formula can be derived for the stable oxygen isotope ¹⁸δ.

Eq. [12] combined with the stable isotope mixing lines [9] and [10] gives

([13]) $$\rm{ ^{14}{a^{SYS}} = 2.05\ {\delta^{18}}O + 97.95\;\;\;for\ model\ P}$$

([14]) $$\rm{ ^{14}{a^{SYS}} = 1.50\ {\delta^{18}}O + 98.50\;\;\;for\ model\ R}$$

where ¹⁴a^SYS represents the recent activity of the estuary, ¹⁴a(est).

The final activity ¹⁴a(fin) for the shell is

([15]) $$\rm{ ^{14}a\left( {fin} \right) {\ =\ ^{14}}a\left( {measured} \right){/^{14}}{a^{SYS}}\left( {est} \right)}$$

The estuary ¹⁴C age can now be calculated as

([16]) $$\rm{ ^{14}C\ age = -8033\ ln\ \left[ {^{14}a\left( {measured} \right) {/^{14}}{a^{SYS}}\left( {est} \right)} \right]}$$

with ¹⁴a^SYS(est) given by [13] or [14] for model P or R, respectively.

b) Observed data

Much information can be derived from a plot of measured δ¹³C vs. δ¹⁸O values for a series of shell carbonates. Figure 3 shows the theoretical mixing lines, i.e. the estuary lines for the stable isotopes. The position of the stable isotope values for the measured shell carbonate determines the mixing ratio between marine- and freshwater reservoirs they were living in, for the two freshwater systems considered: Rhine (R) and other (P).

Figure 3. Selected stable isotope data for shells, showing the estuary lines for the models P and R. Figure 3a shows shells from the Westerscheldt estuary, province of Zeeland; Figure 3b shows shells from Zuid Holland. See text for detailed explanation.

The division between P- and R-water is clearly visible in Figure 3. Shown is ¹³δ vs.¹⁸δ values of a selection of shells subject to estuary mixing of sea- and river water. The selection is that only shells which are not subject to special effects such as atmospheric exchange or uptake of secondary carbonate are shown. These effects will be discussed below (see Figure 4).

Figure 4. Schematic estuary lines (δ¹³C vs. δ¹⁸O) for shells from marine, river and estuary environments and the effect of atmospheric exchange and secondary carbonate.

Figure 3a shows shells from the Westerscheldt estuary, province of Zeeland. It shows measurements of the stable isotopes δ¹³C and δ¹⁸O for historic and recent shells. They are selected from the various datatables. The shells shown are all following the P-line. A few data from the North Sea, just ouside the Zeeland coastline are included. Figure 3b shows a selection of shells from Zuid Holland. In this province the rivers Rhine (R) and Meuse (P) show complex coarse changes during the Holocene (Berendsen and Stouthamer 2002; de Haas et al. 2018).

At present the main flow of the Rhine into the North Sea is near Rotterdam. However, during pre-and protohistoric times the main river flow into the sea was at present Katwijk (near Leiden). This part of the Rhine system still exists as a relatively small river, today called the “Old Rhine.” This water was part of the so-called Limes, the boundary of the Roman empire along the Rhine (van Dinter 2013). The shells shown are older than ca. 6000 years, and younger than ca. 1500 years. These shells all clearly follow the R-line (Rhine), not the P-line (Meuse).

A few datapoints represent unionid shells. This shell is a freshwater taxon, hence they show stable isotope values for riverine water and not estuary water. Indeed, this also shows in Figure 3b; three of the δ¹⁸O values are those of the Rhine (R), indicating that during the periods shown by the dates there was no influence of the sea, or Meuse river (P). These shells include the dated shells GrA-3640 and GrA-3641 (both Maassluis, Zuid Holland, ca. 7200 years old). The third datapoint shown is for a shell from Wijk bij Duurstede near Utrecht (GrN-6326) with a Roman context.

Two other unionid shells are connected with the P-line. One is a recent shell from from the Meuse near Bergen in the province of Limburg (Verdegaal et al. 2005); the other one from Spijkenisse in the Rotterdam region. The shells are not ¹⁴C dated but the Spijkenisse shell originates from an archaeological context dating ca. 4200 years (Versteegh 2009). The stable isotope values represent the river Meuse: δ¹⁸O = –6.5‰ for both. That means that 4200 years ago, the dominant river near Rotterdam must have been the Meuse (P-water).

Figure 3 demonstrates the existence of two clearly distinguished sea/river mixing lines in the Netherlands we call P and R, determining the estuary isotopic values. These estuary lines form the basis of our dating model, together with additional effects determining the final ¹⁴C age of the shell: atmospheric exchange and secondary carbonate.

Figure 4 shows a schematic drawing of the P and R estuary lines, including deviations from the estuary line caused by two different effects: uptake of secondary carbonate from the soil by the shells, and isotopic exchange of CO₂ with the atmosphere. They are indicated by the green and blue lines in Figure 4, respectively. Both effects change the δ¹³C value of the shell, and not the δ¹⁸O value. This makes the additional measurement of the ¹⁸O isotope for shell carbonates crucial, because it determines the position in the δ¹³C/δ¹⁸O plot and thereby also the proper ¹⁴a value, and thus also ¹⁴a^SYS and the ¹⁴C age. This is discussed further below.

The “estuary age” of the shell is –8033 ln(¹⁴a/¹⁴a^SYS). The recent ¹⁴C activity of the estuary ¹⁴a^SYS(est) is calculated by formulas [13] or [14] for model P and R, respectively.

Next a correction term is introduced, representing the difference between the measured δ¹³C value with the value of δ¹³C on the estuary line based on its δ¹⁸O value:

([17]) $${\rm \Delta = {\delta^{13}}{C_{corrected}}-{\delta^{13}}{C_{measured}}}$$

where δ¹³C_corrected is given by the equations [9] and [10] for model P and R, respectively.

The Δ value is used in the correction formulas for atmospheric exchange or secondary carbonate uptake (see Figure 4).

When there is a secondary carbonate effect, the δ¹³C will show below the mixing line because the CO₂ from the soil is depleted in ¹³C. The δ¹³C values originating from the soil organic matter are typically around –25‰. The soil sediment is not recent, hence the ¹⁴C activity of this organic matter is lower than 100%. Therefore, a correction term needs to be subtracted from the measured age.

When there is exchange with the atmosphere, the δ¹³C shows above the mixing line. The ¹⁴C age appears too young, a correction term needs to be added to the measured age. The ¹⁸O in the exchange process does not play a role, because in water the amount of oxygen is much larger than the amount of carbon.

Both effects will be quantified in our model. They are distinguished by their Δ value. The atmospheric exchange correction procedure is used in case Δ < 0; the secondary carbonate correction for Δ > 0. This will be discussed below in sections c) and d).

We use the δ¹⁸O formula for our dating model, because the secondary carbonate effect and the atmospheric exchange effect both change the δ¹³C, and not the δ¹⁸O. Only in a few exceptional cases also the δ¹⁸O is affected; an example is discussed below in section g).

c) Δ < 0: atmospheric exchange correction

The size of the age correction is estimated as 109 Δ. This is determined as follows.

For seawater, the shell carbonate in equilibrium has a recent activity ¹⁴a^SYS = 100% and δ¹³C = 1.5‰. For rivers, the values are ¹⁴a^SYS = 85% and δ¹³C = –10.5‰. The age difference between recent seawater and riverwater is measured as 1300 years, corresponding to a difference of 12‰ in δ¹³C. Thus, 1‰ in δ¹³C corresponds with 109 ¹⁴C years. Thus, the age can be corrected by adding a term 109 Δ years, Δ given by [17].

For the shell age, the atmospheric exchange effect needs to be added to the estuary age given by [16], yielding the final modeled ¹⁴C age:

([18]) $$\rm{ ^{14}C\ age = -8033\ ln\ \left[ {^{14}a\left( {measured} \right) {/^{14}}{a^{SYS}}\left( {est} \right)} \right] + 109\Delta }$$

d) Δ > 0: secondary carbonate correction

In order to quantify the secondary carbonate uptake effect, we introduce a bicarbonate balance equation for the recent activities for the estuary (est) and secondary carbonate (sec):

([19]) $${\rm {(2.5 + F\Delta )^{14}}{a^{SYS}}\left( {tot} \right) = {2.5^{14}}{a^{SYS}}\left( {est} \right) + F{\Delta ^{14}}{a^{SYS}}\left( {sec} \right)}$$

which is based on the following.

The bicarbonate concentration in the ocean is 2.5 mmol/L (Mook 1968, 2006). Extra bicarbonate (the secondary carbonate) is added to the estuary water. This is caused by reactions in the soil which produce HCO₃-. This scales with δ¹³C (in the form of Δ) and an empirical factor F, to be quantified. When Δ = 0, the stable isotopes of the shell are located on the estuary line, then we have simply:

([20]) $$\rm{ ^{14}{a^{SYS}}{\ =\ ^{14}}{a^{SYS}}\left( {est} \right)\;}$$

The factor F has been determined by measurements of carbon properties (chemical and isotopic). For water samples from the Wadden Sea it appears that F = 0.37. For brackish waters from the hinterland, the measurements result in F = 0.88. The latter will be used for stagnant waters, for which we developed a separate “submodel” S which is discussed below. The measurements and interpretations resulting in the above values for F are shown in Appendix B.

All surface reactions (leading to secondary carbonate) are taking place between 0 and 30 cm in the sediment. The sedimentation rate is known to be ca. 2 mm per year, so that the average age will be about 80 years. This corresponds to 1% in the ¹⁴C activity ¹⁴a, so that the recent activity of the sediment can be taken as 99%. For the surface sediment (tidal flat and gullies) this results in

([21]) $$\rm{ ^{14}{a^{SYS}}\left( {sec} \right) = {0.99^{14}}{a^{SYS}}\left( {est} \right)}$$

A special treatment is needed for shell species which do not acquire food by filtering the water, but by grazing organisms (for example Buccinum, Hydrobia/Peringia, and Littorinia species); and also for species applying both filtering and grazing (for example Macoma and Scrobicularia species). See Appendix C for relevant properties of the shell species.

In gullies during low tide algae and seaweed are living in a “recent ¹⁴C” environment (determined by photosynthesis), for which ¹⁴a^SYS = 100%. But the secondary HCO₃ ^– lowers this value to (about) 99%. The biologically active zone is between 0 and approximately 30 cm. Hence here we have mixing of food obtained by filtering with ¹⁴a^SYS = 0.99¹⁴a(est) (R or P) and the surface component ¹⁴a^SYS = 99%, both affected by secondary HCO₃ ^–.

Further, in order to derive the contribution of the secondary carbonate, we introduce a “feeding parameter” γ for the shells. This factor γ ranges between γ = 0 for tidal flat organic matter (shells feeding from the surface only), and γ = 1 for suspended organic matter (shells feeding by filtering, estuary only). This results for the shells in the following properties:

γ = 0 for surface feeding species, γ = 1 for suspension feeding species and γ = 0.5 for species using both (assuming equal proportions).

The recent activity of the secondary carbonate contribution can now be written as a combination of the estuary and tidal flat surface:

([22]) $$\rm{ ^{14}{a^{SYS}}\left( {sec} \right) = 0.99{\gamma ^{14}}{a^{SYS}}\left( {est} \right) + {(1-\gamma )^{14}}{a^{SYS}}\left( {tidal\ flat} \right)}$$

which leads to the following, for model P and R respectively:

model P

([23]) $$\rm{ ^{14}{a^{SYS}}\left( {sec} \right) = 0.99\gamma (1.50{\delta ^{18}}O + 98.50) + 99(1-\gamma ) \approx 99-2\gamma + 2\gamma {\delta ^{18}}O}$$

([24]) $${\rm = {\delta ^{18}}O + 98 (\gamma = 0.5)\;or\ 2{\delta ^{18}}O + 97 (\gamma = 1)}$$

model R

([25]) $$\rm{ ^{14}{a^{SYS}}\left( {sec} \right) = 0.99\gamma (2.05{\delta ^{18}}O + 97.94) + 99(1-\gamma ) \approx 99-1.5\gamma + 1.5\gamma {\delta ^{18}}O\;}$$

([26]) $${\rm = 0.75{\delta ^{18}}O + 98.25 (\gamma = 0.5)\;\;or\ 1.5{\delta ^{18}}O + 97.5 (\gamma = 1)}$$

Some shells are carnivorous, for these γ = 0. This leads for both models P and R to

([27]) $$\rm {^{14}{a^{SYS}}\left( {sec} \right) = 0.99}$$

All parts of equation [17] are now determined, which enables calculation of the total secondary carbonate contribution:

([28]) $$\rm {^{14}{a^{SYS}}\left( {tot} \right) = [{2.5^{14}}{a^{SYS}}\left( {est} \right) + F{\Delta ^{14}}{a^{SYS}}\left( {sec} \right)]/(2.5 + F\Delta )}$$

Next we can calculate the final activity for the dated sample:

([29]) $$\rm {^{14}{a^{SYS}}\left( {fin} \right) = 100 \left[ {^{14}a{/^{14}}{a^{SYS}}\left( {tot} \right)} \right]}$$

where ¹⁴a is the measured activity. Finally this leads to the modeled age:

([30]) $$\rm {^{14}C\ age = - 8033\,ln{[^{14}}a\left( {fin} \right)/100]}$$

where ¹⁴a^SYS(fin) is given in %.

e) Stagnant waters: model S

Shells living in stagnant seawater (only Cardium glaucum species can survive in this water) form a special category for the treatise of secondary carbonates.

The ¹⁸δ value can not be used as a tracer for stagnant waters. It is influenced by evaporation and precipitation as a function of time. This also means we can not use the same Δ as before, since ¹³δ_corrected depends on the ¹⁸δ value.

For the stagnant water case, we use an alternative Δ here called Δ* (* to avoid confusion):

([31]) $${\rm \Delta^{*}{\ =\ ^{13}}{\delta _{sea}}{-^{13}}{\delta _{measured}} = 1.5{-^{13}}{\delta _{measured}}}$$

since the ¹³δ value for the sea is 1.5‰ (see Table 1).

The situation discussed here is the mixing between seawater and the secondary carbonate. The case discussed above in section c) dealt with mixing of estuary water and the secondary carbonate. Taken together, the formula now translates into

([32]) $$\rm {^{14}{a^{SYS}}\left( {tot} \right) = [{2.5\ ^{14}}{a^{SYS}}\left( {sea} \right) + 0.88\ \Delta^{*}\ {^{14}}{a^{SYS}}\left( {sec} \right)]/(2.5 + 0.88\Delta^{*})}$$$${\rm = [250 + 0.88\Delta ^{*}\ {^{14}}{a^{SYS}}\left( {sec} \right)]/(2.5 + 0.88\Delta^{*})}$$

with ¹⁴a^SYS = 100% for the North Sea.

The secondary carbonate contribution is quantified as ¹⁴a^SYS(sec) = 0.97 = 97%. This value is based on well dated stratigraphies of peat dates, so-called vegetation horizons and marine clay deposits from the province of Groningen. The marine clay contains shells deposited during storm surges which lived in the remaining stagnant pools. The data and a full analysis leading to the number 0.97 = 97% is shown in Appendix D.

For model S, the formulas [29] and [30] derived above can be applied to calculate the final modeled ¹⁴C age.

f) Summary of formulas

The essential formulas for the various (sub)models are summarized in the textboxes below. The activities ¹⁴a are in %, and the δ values in ‰. Examples of worked calculations are given in the Appendix G.

model P

model R

model S

g) Special mixing situations

The model is expanded in order to deal with complex mixing situations, often changing during (pre/proto)historical times. Most prominent examples are the northern Rhine branch known as Oer-IJ in North Holland (Vos 2015) and the Rhine/Meuse delta in Zuid Holland (Berendsen & Stouthamer, 2002; Pierik et al. 2018). Detailed geophysical analysis using our model will be published in a more applied local journal.

This leads to (sub)models called “L” (lake) and “Z” (Zuiderzee). The model calculations are shown for “P” and “R” in the main tables together with the administrative information; the final dates (modeled and calibrated) are given for “L” and “Z.” Detailed calculations for “L” and “Z” are given in separate tables. These special submodels only apply to a limited amount of shell dates.

Applying the model: illustrative examples

The applications of our model for the physical geography and archaeology of the Netherlands are many. In this paragraph selected examples of model applications will be discussed. These are illustrative and show the validity of the model. The examples concern estuary shells (a), shells subject to atmospheric CO₂ exchange (b) or secondary carbonate uptake (c), and others (d–f). All dates in this section are rounded to 5.

a) Estuary shells

The oldest Holocene shells within the recent coastline are dated about 7000 ¹⁴C years old. These are Cardium shells from the locations Beverwijk (N. Holland), Maasdijk (Z. Holland) and Kerkwerve (Zeeland); see Table 2, which also contains dated associated and contemporaneous terrestrial samples.

Table 2. Examples of modeled estuary shell dates compared with contemporaneous samples

The locations represent Early Holocene estuaries from three different locations. Beverwijk is located west of Amsterdam at the present coastline; at the time of the shell, it was part of a large estuary of the Vecht river (see Figure 2). The Beverwijk shell is from a depth of 16.98 m NAP (New Amsterdams Peil, the Dutch Ordnance Datum). It is deposited on top of a peat layer known as base peat, which is is dated as 7260 ± 60 BP (GrN-15928). The uncorrected age (i.e. not normalized, not modeled) and conventional age (after subtracting 400 years) appear way too old. Applying our model P yields a date of 7155 years, which is consistent with the peat date: it is younger. The Maasdijk shell originates from a depth of 26.2 m NAP, in a deep marine gully. The location is between Rotterdam and the coast. It should date younger than the base peat deposits in the region; see the cross section in Vos (2015, 242). Peat stratigraphy from nearby Rotterdam has been dated, the relevant date is 7160 ± 45 BP (GrA-32101). See Cohen and Hijma (2008) for a treatise of the Rotterdam (Blijdorp) excavation, and Hijma and Cohen (2011) for a palaeogeographic description of the region.

The uncorrected shell age, 7585 years is too old, as is the normalized age with additional subtraction of 400 years. Applying our model R yields a date of 6880 years, which is consistent with the geography. Kerkwerve is located on the island of Schouwen-Duiveland in the province of Zeeland. The shell represents the old Scheldt river course. Organic samples from the region date 7280 ± 45 BP (Anna Jacoba polder on Tholen, peat, GrN-16936) and 7185 ± 40 BP (Korendijk, charred seeds, GrA-49578). Stratigraphically, the shell should date younger. Also here, the not corrected date is significantly older; the same is true for the normalized age in BP, including subtraction of 400 years. The P-modeled age is 6835 years which does fit the stratigraphy.

The value of 400 years for the reservoir effect in the examples above is obviously not correct, that applies only to a 100% marine shell habitat. We emphasize this here because this value has probably been used in the users community. These examples (as well as the ones following below) illustrate the first successful method to quantify the reservoir effect that fits the stratigraphy, based on geophysical observations and ¹⁴C dates of terrestrial samples.

b) Leijduin: shells subject to the atmospheric exchange effect

Leijduin is located in the province of North Holland. Measurements for Cardium coastal shells from this location yield the isotopic ratio values δ¹⁸O = –4.53‰ and δ¹³C = –2.58‰, and ¹⁴a = 58.80% (GrN-14203). The uncorrected ¹⁴C age is –8033ln(¹⁴a) = 4265 years.

The estuary age is calculated using ¹⁴a^SYS(est) = 2.05 ¹⁸δ+97.94 = 88.65% which results in an estuary age of 3298 years. The shell lived in an environment with no uptake of secondary carbonate from the soil. The stable isotope datapoint is located far above the estuary mixing line. The estuary line ¹³δ = 1.64 ¹⁸δ-0.14 for the P model predicts δ¹³C = –7.57‰, revealing a large difference of Δ = 4.99‰ caused by isotopic exchange with the atmosphere. This results in a shift in age of 109Δ = 543 years, yielding a modeled age of 3298+543 = 3841 years.

The δ¹³C value is strongly enriched. Other large enrichents are observed for GrN-7774 and 7775 (Alkmaar), GrN-9042 (Heiloo), GrN-10169 (Nijenburg) and GrN-10170 (Kooibrug).

The Leijduin shell represents the most extreme case in our dataset. The estuary effect is close to 1000 years in one direction, and the exchange effect more than 500 years in the other direction. Observe that for a shell not subject to uptake of secondary carbonate, the equation ¹⁴a^SYS = ¹⁴a(est) = 1.25¹³δ + 98.25 = 94.90%, resulting in a ¹⁴C age of –8033ln(58.80/94.90) = 3845 years, which is practically the same as the number 3841 above. It also shows a reservoir effect of 4265–3845 = 420 years. All together, this analysis shows the validity of our correction model for ¹⁴C dates.

c) Shells subject to secondary carbonate and the sea level

Figure 5 shows ¹⁴C dates of shells from the province of Zeeland, plotted as a function of depth. They can be compared with the sea level curve of Meijles et al. (2018), known as the MSL curve. The latter is obtained by ¹⁴C dating of peat deposits and is valid for the part of the North Sea discussed here. The peat record does not extend beyond a depth of 12m.

Figure 5. Age/depth plot of shells from Zeeland, shown together with the sea level curve of Meijles et al. (2018) which is based on peat dates. The Scrobicularia shells follow the sea level curve (see text for details).

Figure 5 shows uncorrected (measured) and corrected (modeled) ¹⁴C dates, using model P for the ¹⁴C dates and IntCal20 for calibrated dates. The model makes the shells younger. The black open circles show the uncorrected dates; the filled circles the modeled Scrobicularia (scp) shell dates, and diamond symbols for selected Cardium (car) shell dates.

For the deepest part (below 14 m) of the shell curve, only a few measurements for car shells (no scp shells) are available. Observe that the sea level curve based on peat (the MSL curve) represents dated plant remains, whereas dates of shells can provide a more direct measurement of the sea level at low tide. The MSL peat curve is shown as a green dashed line (see Figure 10 in Meijles et al. 2018). The sea level curve for shells (based on model P) is shown as a purple dashed line drawn through the datapoints and is located below the peat curve. This analysis shows that from our shell database, sea levels can be reconstructed. Our model is shown here to provide a valid explanation for the measured ¹⁴C dates.

The shell record contains some shells from much larger depth. Also, the shell dates as a function of depth represent the sea level, albeit of a different nature than the MSL curve. It concerns only shells that live at or near the dry surface, about 0–30 cm depth during low tides. This concerns Scrobicularia (scp) shells, and often Cardium edule (ced) shells. The latter can also be found in deep waters, as is also the case for Macoma (mac) shells. The near surface shells should show below the MSL curve; they can not be located above this curve. The dates of the scp shells do fit this criterium, after correction as estuary shells including the secondary carbonate effect.

An interesting observation is the scp shell that shows well above MSL curve (GrA-3528). As it happens its ¹⁴C date (modeled age 7285 ± 130 ¹⁴C years; the calibrated age 8375–7860 calBP) corresponds to the well known Storegga tsunami event in the North Atlantic, including the North Sea (Walker et al. 2020 and references therein). The shell is not in situ and must have been redeposited later.

d) Monster: an example of stagnant water

Monster is a small town, located south of Den Haag (The Hague) at the coastline. The natural landscape is a clay deposit on top of peat; the peat layer (called “Hollandveen”) throughout the region has known ¹⁴C dates of about 2600–4200 BP from top to bottom. Below this peat an older layer of clay was deposited.

The Monster site is an archaeological excavation, showing Iron Age activity; peat was reclaimed in the past by digging holes. One of these holes contained a wooden artifact. Later, this hole was filled with marine deposits and became brackish, including Cardium glaucum (cgl) shells that favour such environments. Again later, the pit was covered by phragmites-sedge peat (“rietzeggeveen;” de Jong, ¹⁴C sample submission form RGD; van der Valk [1988]).

Hence, the stratigraphy (from top to bottom) is peat (rietzeggeveen), shells, artifact, peat (Hollandveen) and clay. The bottom clay layer also contained marine shells. The artifact, shells and top peat deposit have been ¹⁴C dated. For a schematic drawing of site and sample locations see Figure 6.

Figure 6. Schematic drawing of an archaeological site at Monster (Zuid Holland), illustrating the validity of our shell dating model for stagnant waters. The peat date GrN-15418 is 2300±20 BP, the wood date GrA-882 is 2360±60 BP. The modeled shell ages are 2310±60 (GrN-15416), 2540±65 (GrN-15417), 4360±50 (GrN-15108) and 4365±60 (GrN-15109).

The phragmites-sedge peat dates to 2300±20 BP (GrN-15418), and the wooden artifact to 2360±60 BP (GrA-882). Both dates correspond to ca. 350 BC. Two Cardium glaucum (cgl) shells from the marine sediment were dated; they must be younger than the wood artifact. But both GrN-15416 (uncorrected age 2495±60) and 15417 (uncorrected age 2720±65) date too old when normalized, and a reservoir effect of 400 years subtracted. After correction using model S the ¹⁴C age of the shell GrN-15416 is 2310±60, confirming the expected age; GrN-15417 dates 2540±65 and is still too old. The latter observation can possibly be explained as follows. Clay sediment is partly mixed with peat particles affecting ¹⁴a(sec) with ¹⁴a(sec) =97%, causing an “old peat effect”. Based on the average age of the Hollandpeat and assuming a contaminating amount of 10%, we obtain ¹⁴a(sec) =94% and the date becomes 2355 modeled ¹⁴C years.

This example shows the validity of our model S for stagnant waters. It is important to note that Cardium glaucum shells where often not submitted for dating by the RGD (Geological Survey of the Netherlands), because “they always appear too old” (W.H. Zagwijn, unofficial statement on submission form), either in BP and subsequent marine reservoir correction, or not normalized for fractionation. Our model now provides an explanation for this interesting observation made decades ago.

e) The province of Groningen

A particular observation is that throughout the northern part of the Netherlands, the then existing peat deposit was covered with marine clay during the 9th century BC. The marine deposit is interpreted as the result of several storm surge events during 840–760 BC. The peat deposits are extensively dated by ¹⁴C; throughout the northern Netherlands, the typical ¹⁴C date for the toplayer of this base peat is 2600 BP (see also Appendix A). The marine deposit contains shells, many of which are dated by ¹⁴C. We applied our model to derive calender dates, focussing on the province of Groningen from which we have collected and analyzed dozens of samples.

An overview of shells from the province of Groningen is shown in Figure 7a. The stable isotope ratios δ¹³C and δ¹⁸O are shown, with the estuary lines for model P and R (see also Figure 4). Figure 7a shows the stable isotopes of all shells (δ¹³C vs. δ¹⁸O), categorized by species Cardium edule (ced), Cardium glaucum (cgl), Macoma (mac), Scrobicularia (scp) and others. The same data are shown in Figure 7b, but now categorized by their Δ value which is the difference between measured and estuary δ¹³C value used in our model. The plots represent shells from all (Holocene, but not recent) ¹⁴C ages of the province.

Figure 7. Shell stable isotope data for the province of Groningen. The dashed lines show the estuary mixing lines for precipitation fed rivers (P) and the Rhine (R). a) shells categorized to shell species: Cardium edule, Cardium glaucum, Macoma and Scrobicularia; b) shells categorized to their Δ values (see text).

Considering uncertainties involved, we have chosen Δ<–0.25‰ for shells subjected to the atmospheric exchange effect, Δ>1.0 ‰ for shells subjected to the secondary carbonate effect, and Δ values Δ ≈ 0 ‰ (between –0.25‰ and +1.0‰) for shells only subjected to estuary mixing.

The dataset for the province of Groningen contains one “outlier,” a freshwater species with δ¹³C =–4.42‰ and δ¹⁸O =–3.68‰ which is not part of our modeling. It is shown for completeness. All others are marine shells.

As can be seen in Figure 7b, the data with Δ>1 concern only Scrobicularia, Macoma and Cardium glaucum shell species. The data with Δ≈0 are Cardium shells, positioned on or close to the precipitation model line. There are only a few datapoints with Δ<1. It is obvious thet the datapoints with Δ>1 are dominating; the region was characterized by tidal flats and stagnant waters, resulting in (often large) secondary carbonate effects.

The Groningen coastline was, and still is characterized by an extended tidal flat. Today the land is protected by a high seadike. In historical times, land was continuously reclaimed and protected by dikes still visible in the landscape. During prehistory, people lived on “terpen” (mounds). The land was regularly flooded by storm surges, often catastrophic. The flooding events left their traces, as visible in shell middens with especially Cardium (car, ced, cgl), Scrobicularia (scp) and other species, but not Cardium glaucum in stagnant water, cgl(S). Macoma (mac) and Scrobicularia (scp) are shells living in gullies in the tidal flat, which are rich in organic matter. Hence, they show large secondary carbonate effects, as is clear from Figure 7.

Also, cgl(S) shells are subject to even larger secondary carbonate effects, because seawater is flooding organic rich vegetations. Isotopically, the three species are very similar. However, the mac and scp species require running water conditions; they can not live in stagnant water, whereas for the cgl(S) shells this is the favorite habitat. During storm surges seawater flooded the hinterland and became stagnant. When shell brood of Cardium shells (car, ced) is involved, they will develop into cgl(S) in these circumstances. They die when water becomes too fresh or disappears by evaporation, in fact the living period of the shell is an indication for an “event.” Their isotopic behavior is described by our model S for stagnant waters. See also Appendix C for information on the various shell species discussed in this article.

A storm surge with resulting stagnant water as described above occurs during a relatively short time. Therefore, the cgl(S) shells usually are an indication for an “event.” Note that there are also other types of storm surges, dumping shell middens on the tidal level, even above sea level. These contain various species, but not cgl(S); these are redeposited from the sea bottom.

Following the retreat of the brackish seawater, pastures developed recognized by so-called vegetation horizons. These are characteristic thin dark layers of organic material in Late Holocene clay sediments. Such horizons represent moist meadows (Taayke et al. 2022). These horizons provide an important time marker in the region, for instance the Early Iron Age (ca. 2500 BP). Younger horizons are also observed.

The storm surge period during 840–760 BC coincides with a time characterized by significant climate change, rising groundwater tables and drowning of land, and migration of inhabitants known in the paleobotanical and archaeological record (van Geel et al. 1996, 1998; van Zijverden 2017). It also coincides with a large increase in the Δ¹⁴C record, part of the millennial scale oscillation in solar cosmic ray production and apparently also in the North Atlantic climate (Bond et al. 2001). The 840 BC event also shows in some shells from the coast of Noord Holland and Zuid Holland; in Zeeland it is not observed—here, no shells for this age are found or measured.

As a summarizing example of model S, Table 3 shows a direct comparison of shell dates and those of associated peat. The peat samples are representative for a large series of peat dated througout the region.

Table 3. Examples of shell dates (using model S) compared with associated peat dates

The shells are deposited on the peat and should therefore be younger. The modeled dates are younger; the not-modeled (uncorrected) dates are all too old. Hence, the model is in accordance with the stratigraphy, and the not-modeled dates are not.

See also Appendix D, which contains more information on the peat deposits leading to the value 97% for the recent value of the secondary carbonate ¹⁴a^SYS(sec), which is used for model S.

f) Complex mixing situations in the Zuiderzee region

The shells from Zuiderzee (Lake IJssel, Flevoland; see Figure 1) represent complex mixing situations, depending on time period and location. The same is true for the so-called Oer-IJ region (Vos et al. 2015; Weerts et al. 2002).

Our pre-modern shell dataset for the province of Flevoland (the Lake IJssel region) is limited in size—only 9 shells have been dated. Nevertheless, their isotope values illustrate the model in terms of past geophysical developments. Figure 8 shows the measured isotope rations δ¹³C and δ¹⁸O for two distinguished groups, according to their age.

Figure 8. Shell stable isotope data for the IJsselmeer (Lake IJssel) region. The dashed lines are the estuary mixing lines for precipitation fed rivers (P) and the Rhine (R). Two groups can be distinguished: “old” dates (red) and young dates (blue). This can be explained by the breakthrough of Rhine river water (via the IJssel) to the sea during the first millennium AD.

For a palaeogeographical description of the IJsselmeer region, we refer to Vos and de Vries (2016), Grondboor & Hamer (2016), van den Biggelaar (2017) and van den Biggelaar et al. (2015). The Rhine river became connected to the Lake IJssel during the first millennium AD. When exactly is not known, Cohen et al. (2009) conclude around 550 AD, and Vos (2015) states ca. 800 AD.

The 5 older dates (ca. 3500 years, thus much older than the Rhine breakthrough) are shown in red. They do not show influence of Rhine water, as it should be for shells of this age. The 5 shells appear to belong to 2 groups; 2 (GrN-11498, 19513; C in Figure 8) are on the precipitation line, and 3 (GrN-7082,10623,18358; D in Figure 8) are shifted to more positive δ¹⁸O values. The location of shells “C” is close to the mouth of the Eem river into the lake, which is probably the reason they are not subject to secondary carbonate. The shells “D” are below the estuary line, caused by the uptake of secondary carbonate.

The Rhine component should be only visible in the younger set of dates. They are shown in blue in Figure 8. They are from two locations (Hindeloopen and Medemblik), with two different shells analysed for each location. The two with the most positive δ¹³C values (A in Figure 8) are Cardium glaucum (cgl) shells. They show above the Rhine line, subject to exchange with the atmosphere, and their isotope values are consistent with those of modern shells. The two with the more negative δ¹³C values (B in Figure 8) are Scrobicularia (scp) shells. The values are below the Rhine line, so these are clearly subject to the secondary carbonate effect.

For the younger group we developed model Z, with specific calculations for individual shells. The modeled dates are shown in the datatables. Here we show as an example a detailed calculation for two shells from Hindeloopen, a small town located at the Zuiderzee coast of Friesland.

For the Zuiderzee, Eisma, Mook and Das (1976) analyzed the stable isotopes (¹³C and ¹⁸O) for historic shells obtained from museum collections. They were collected before 1932, the year that the Afsluitdijk was completed, which transformed part of the sea into a lake. For these shells no context is known, except for its general location along the coastline. A limited amount was ¹⁴C dated later by the conventional method.

Below, we use the subcript “m” for marine water, “f” for freshwater and “s” for shell carbonate. Eisma et al. (1976) determined the Chlorinity for Cardium shells; the [Cl] value determines the number of growth bands R (known as “shell ribs”). The relation between both quantities is given by the equation [Cl]_s = 2R-31‰ (curve “E” in Figure 1 in Eisma et al. [1976]).

The mixing line (¹⁸δ_m-¹⁸δ_s)/([Cl]_m-[Cl]_s) between Zuiderzee and North Sea leads for the measured properties to the relation between the ¹⁸δ value and chlorinity of the shells:

([33]) $${\rm {\left[ {Cl} \right]_s} = {2.75^{18}}{\delta _s} + 16.52}$$

The measured isotope values and Chlorinity calculated using [33] for the Hindeloopen shells are shown in Table 4. The marine percentages are calculated by using [Cl]_s/19.3.

Table 4. The measured data for the isotopes ¹³C, ¹⁸O and ¹⁴C and Chlorinity data for two shell species from the Zuiderzee

The Cardium (cgl) shell is subject to the atmospheric exchange effect; the Scrobicularia (scp) shell, is subject to secondary uptake of older carbon from the soil.

The mixing equation for marine- (relative amount x) and freshwater (amount 1–x) is

([34]) $${\rm \left( x \right){(^{13}}{\delta _m}) + (1-x){(^{13}}{\delta _f}){ = ^{13}}{\delta _s}}$$

The marine component is 38% =0.38, hence for ¹³C, δ_f = (δ_s –1.5x)/(1–x) which gives ¹³δ_f = –3.95‰; for ¹⁸O, δ_f = (δ_s –x)/(1–x) resulting in ¹⁸δ_f = –6.02‰ for GrN-12024.

We can now derive the initial ¹⁴C concentration ¹⁴a^SYS for the freshwater component for the shell, using the equation [12]: ¹⁴a_f = 1.25¹³δ_f + 98.125 = 93.19%. The isotope parameters for the estuary water can be determined from the mixing equations for seawater/freshwater [11] resulting in ¹⁴a(est) = 0.97¹⁸δ_s + 99.03 and ¹³δ(est) = 0.78¹⁸δ_s + 0.72.

For the scp shell, we can determine now also the secondary carbonate contribution ¹⁴a(sec), using mixing relations between estuary and tidal flat (wad), for which ¹⁴a^SYS = 100%. Further we assume equal contributions from wad and estuary (1:1 mixing).

Next, we need to establish the relationship between ¹⁴a(sec) and ¹⁸δ_s, as follows.

For ¹⁸δ_s = –3‰, ¹⁴a(est) = 96.12% yielding ¹⁴a(sec) = [(96.12+100)/2](0.99) = 97.08%;

for ¹⁸δ_s = 1‰, ¹⁴a(est) = 100% resulting in ¹⁴a(sec) = [(100+100)/2](0.99) = 99%.

The mixing line between these 2 points (¹⁸δ_s, ¹⁴a(sec)) = (–3, 97.08) and (1, 99) then finally gives for the secondary effect ¹⁴a(sec) for this specific shell: ¹⁴a(sec) = 0.48¹⁸δ_s + 96.60

Finally, the modeled ¹⁴C dates are calculated, using equations [12] and [28] for GrN-12024 and GrN-12025, respectively:

¹⁴a^SYS = 1.25¹³δ_s + 98.125 = 95.775% and ¹⁴C age = −8033 ln[¹⁴a_s/¹⁴a(sys)] = 440 years;

¹⁴a^SYS = [2.5¹⁴a(est) + 0.37Δ¹⁴a(sec)]/[2.5 + 0.37Δ] = 95.90% and ¹⁴C age = –8033 ln(¹⁴a_s/¹⁴a^SYS) = 603 years.

Estimation of model uncertainties

a) General remarks

The ¹⁴C dates are the result of a physical measurement, for which uncertainties (error bars) are determined. For a Gaussian probability distribution, these are given as sigma σ values in BP ± σ, representing 68.3% probability. The calibration curve, the relation between ¹⁴C dates and the calender, is not a linear function. Therefore, the calibrated dates are not Gaussian, and here the probability distribution is given as a 68.3% dating range (corresponding to 1σ in BP) or a 95.4% dating range (corresponding to 2σ).

The uncertainties introduced by handling and interpreting the measured ¹⁴C dates for aquatic samples are not straightforward to quantify. This concerns variations in the size of the reservoir effect, and uncertainties in the applied model parameters.

For this reason, the calibrated dates are reported as 95.4% ranges, taking the complete range and ignoring the effect of multiple peaks in the distribution, when present. This is a conservative and pragmatic approach. It is clear that in general the model works—it provides an explanation of ¹⁴C dates that appear wrong, when the model is not applied. Prime examples are given in the chapter showing applications of the model. This holds for practically all our dates, which represents an unprecedented number of dated shells.

Here we argue that the uncertainties in the values of the main model parameters used do not significantly change the outcome of the calibrated dates. The resulting increase of larger uncertainties for the 95.4% probability range in calender years is limited in almost all cases. Below follow estimates for the size of the uncertainties, calculated by varying the values of the most significant parameters used in the model.

b) The recent activity of the freshwater reservoir

The main parameter is the recent activity of the riverwater, which is determined to have the value ¹⁴a^SYS = 85%. This is based on measurements and other observations, discussed in Appendix A. There we argue that a reasonable assessment for the uncertainty in the ¹⁴a^SYS value for river water is about 1%.

The influence of this uncertainty on the modeled age can be best shown by performing calculations using 84% instead of the model value 85% in equation [11]. For a shell with δ¹⁸O = –2‰, the calculated result for the ¹⁴C activity using model P is then ¹⁴a = 93.4‰. For the chosen model parameter value ¹⁴a_f = 85%, the calculated activity will be ¹⁴a = 93.8‰. Thus, the difference between the two calculations is 0.4% in ¹⁴a. Hence, a deviation of 1% in the number 85% is equivalent to about 30 ¹⁴C years.

This effect of the ¹⁴a_f value deviating from 85% depends mainly on the δ¹⁸O value (or alternatively Chlorinity) of the shell. Calculations show that the effect becomes larger for more negative δ¹⁸O values. The result of 30 years as calculated above can be considered as a typical value for our dataset.

In a similar way the effect of variations in the input numbers for ¹⁸δ_f can be estimated. These are assessed as –9.0 ± 0.2% and –6.3 ± 0.2% for model R and P, respectively. In this case, the effect for the modeled ages is not larger than 20 ¹⁴C years.

Summarizing, the value 85 ± 1% for the recent activity (and its uncertainty) of the freshwater is shown to be valid for the Rhine/Meuse delta for the part of the Holocene investigated. It is not necessarily a valid number for other regions. See Appendix A for more details.

c) The secondary carbonate

The factor F is used in the formula to calculate the secondary carbonate contribution. It shows two different values for two different environments. These values are determined by a combination of experiments and analysis of data, in particular ¹⁴C dates of peat (see Appendix B for details). Including uncertainties these values are assessed as F = 0.37 ± 0.10 for outside the dikes, and F = 0.88 ± 0.10 for inside the dike. The effect on the calculated ¹⁴C ages is negligible (<10 years).

d) Illustrative example

The assessment of the uncertainty in the modeled age is best illustrated using an example. The shell GrM-15351 is considered a high precision date, with a standard deviation of the measurement of 15 BP. The modeled age is calculated as 2165 ± 15 ¹⁴C years. For this shell, the calibrated probability distribution of the age is bimodal; the calibration graph graph shows 2 peaks, around 2280 and 2140 calBP. There is no probability in the middle (the median), around 2200 calBP. Both peaks taken together results in a calibrated age range 2295–2120 calBP (68.3%), or 2300–2100 calBP (95.4%). We use the 95.4% probability, and final dates are rounded to 5.

This high precision date results in a modeled calibrated age range of about 2 centuries. When we increase the measurement uncertainty from 15 to 20, the resulting modeled and calibrated date range will not be significantly different. This example can be considered an extreme case, but it does illustrate the relative insignificance of the uncertainties in the model parameters.

e) Final notes

We note that uncertainties also apply to the “standard model” using a 400 year reservoir effect. The 400 year value is not a holy number either, also here uncertainties can be larger than the standard deviation of a high precision ¹⁴C date.

The calibrated numbers are taken together as one range, also when it includes wiggles and/or multiple solutions. In most cases, these are not very relevant for aquatic samples anyway, since these are “damped” compared with terrestrial samples. The latter reflect changes in the atmospheric ¹⁴C content (moreless) instantaneously.

Adding all uncertainty contributions σ_i quadratically following the standard statistical formula σ_tot = √[1/Σ(1/σ_i ²)] yields a final total uncertainty in the modeled age not larger than 2σ, this σ being the reported measurement error. The 95.4% callibrated age probabiity distribution we therefore consider as a reasonable pragmatic approach for the error in the modeled age. This is supported by the modeled dates which fit the stratigraphy taking into account the uncertainties, whereas not modeled dates mostly do not. Uncertainties remain, but based on our large dataset we can derive parameters that explain the observed modeled dates and their uncertainties with confidence.

Conclusion

A model is presented to assess reservoir effects for shells fom the North Sea, rivers and complex mixing waters such as estuaries and lakes during the Holocene in the Netherlands. The model and background information leading to it are described and discussed, and all data (measurements and calculations) are included.

The shells are all dated in Groningen during the last 7 decades. A grand total of 1116 shells are dated and presented in datelists. The whole covers unique aspects: the size of the dataset (many hundreds of shells), representing a large density for a relatively small region which is well studied in terms of palaeogeography and archaeology. For most shell carbonates, the three isotopes ¹⁴C, ¹³C and ¹⁸O are measured, thus including ¹⁸O which is/was not standard practice, depending on period of measurement and/or laboratory. The two stable isotopes combined are tracers for the past environment the shells lived in, which is crucial information to derive the size of the reservoir effect.

For the Netherlands, the size of the reservoir effect is difficult to assess, as the shells often lived in an environment of mixed marine- and river waters. River water occurs in two main categories, distinguished by ¹⁸O: the Rhine which is dominant, and other rivers. This leads to two estuary mixing lines between the North Sea and rivers. The stable isotopes of the shell carbonate are also indicative for additional processes, such as uptake of secondary carbonate from the soil by shells, and exchange of C isotopes between atmosphere and water. Extensions of the main model deal with special cases such as pools of stagnant water and lakes.

The model leads to an assessment of the recent ¹⁴C activities of the system the shells lived in, which we call ¹⁴a^SYS. The measured ¹⁴C activities relative to these ¹⁴a^SYS values determine the ¹⁴C age of the shells and include the reservoir effect. This way we circumvent normalizing to δ¹³C = −25‰, i.e. the terrestrial timescale and the subsequent correction for reservoir effects.

We emphasize that the modeled dates are not in BP, because that would include normalization to δ¹³C = –25‰, the conventional correction for isotope fractionation. Perhaps superfluous, it is noted that the other elements of the convention (the oxalic standard and Libby half-life) are used. Our modeled dates can be calibrated using the present calibration curve IntCal20.

The model is based on the three isotopes ¹³C, ¹⁴C and ¹⁸O for shell carbonate. The model uses ¹⁴a^SYS = 0.85 as the recent value for rivers. The ¹³C isotope ratio for riverine shell carbonate is δ¹³C = −10.5‰. For ¹⁸O, the Rhine shows δ¹⁸O = −9.0‰, for the other rivers δ¹⁸O = –6.3‰ applies. For the North Sea, the values are δ¹³C = 1.5‰ and δ¹⁸O = 1.0‰.

We recommend that the ¹⁸O isotope should receive more attention by the user community. The position of the ¹⁸O isotope in the δ¹³C/δ¹⁸O plot determines primarily the ¹⁴C value of a mixed reservoir, and thus the proper shell age. The present article concerns the Rhine/Meuse delta in the Netherlands. We recommend other delta regions to be investigated using the same or similar methods, using the proper regional values for the three isotopes.

We also recommend that in dating reports it is completely clear which ¹⁴C property is reported: ¹⁴a (i.e. not normalized for fractionation) or ¹⁴a_N (i.e. normalized, equivalent to F¹⁴C). For ¹⁴C dates, it should be clear wether it is calculated from ¹⁴a or ¹⁴a_N. Formally, only the latter leads to dates in BP. The first includes reservoir effects, the latter requires correction for reservoir effects.

For the Netherlands, the validity of the model is shown by many examples. Only a limited selection of these is discussed here; others will follow in more local publications. The modeled dates are consistent with stratigraphy; not modeled dates are inconsistent with stratigraphy.