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 Reference Berendsen and Stouthamer2002; de Haas et al. Reference de Haas, Pierik, vsn der Spek, Cohen, van Maanen and Kleinhaas2018).

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 Reference van Dinter2013). 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. Reference Verdegaal, Troelstra, Beets and Vonhof2005); 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 Reference Versteegh2009). 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 Reference Mook1968, Reference Mook2006). 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 Reference Vos2015) and the Rhine/Meuse delta in Zuid Holland (Berendsen & Stouthamer, Reference Berendsen and Stouthamer2002; Pierik et al. Reference Pierik, Stouthamer, Schuring and Cohen2018). 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.

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 (Reference Vos2015, 242). Peat stratigraphy from nearby Rotterdam has been dated, the relevant date is 7160 ± 45 BP (GrA-32101). See Cohen and Hijma (Reference Cohen and Hijma2008) for a treatise of the Rotterdam (Blijdorp) excavation, and Hijma and Cohen (Reference Hijma and Cohen2011) 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. (Reference Meijles, Kiden, Streurman, van der Plicht, Vos and Gehrels2018), 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. (Reference Meijles, Kiden, Streurman, van der Plicht, Vos and Gehrels2018) 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. Reference Meijles, Kiden, Streurman, van der Plicht, Vos and Gehrels2018). 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. Reference Walker, Gaffney, Fitch, Muru, Fraser, Bates and Bates2020 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 [Reference van der Valk1988]).

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. Reference Taayke, Streurman and van der Plicht2022). 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. Reference van Geel, Buurman and Waterbolk1996, Reference van Geel, van der Plicht, Kilian, Klaver, Kouwenberg, Reynoud-Farrera and Waterbolk1998; van Zijverden Reference van Zijverden2017). 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. Reference Bond, Kromer, Beer, Muscheler, Evans, Showers, Hoffmann, Lotti-Bond, Hajdas and Bonani2001). 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. Reference Vos, de Koning and van Eerden2015; Weerts et al. Reference Weerts, Cleveringa and Gouw2002).

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 (Reference Vos and de Vries2016), Grondboor & Hamer (2016), van den Biggelaar (Reference van den Biggelaar2017) and van den Biggelaar et al. (Reference van den Biggelaar, Kluiving, Bohncke, van Balen, Kasse, Prins and Kolen2015). The Rhine river became connected to the Lake IJssel during the first millennium AD. When exactly is not known, Cohen et al. (Reference Cohen, Stouthamer, Hoek, Berendsen and Kempen2009) conclude around 550 AD, and Vos (Reference Vos2015) 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 (Reference Eisma, Mook and Das1976) 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. (Reference Eisma, Mook and Das1976) 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. [Reference Eisma, Mook and Das1976]).

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.

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.