Assessment of the Performance of 210Pb-Based Dating Models with a Challenging Sediment History in Maryport Harbour (UK)

Abstract

The ²¹⁰Pb-based method is used for absolute age determination in recent sediments (<150 years). Different assumptions are possible, leading to different models and chronologies. The evaluation of the capacities and limitations of the models in challenging sedimentary scenarios is of broad interest to the scientific community, and this is the aim of the present work. The performance of the classical models, CFCS and CRS, and the novel TERESA is assessed with a 2 m long, high-resolution core sampled in Maryport Harbour, UK, by using raw data from the literature. It was affected by dredging, shortening, and by the anthropogenic impacts of radionuclides released by the Sellafield nuclear plant and the phosphate industry in Whitehaven and was considered non-datable by the ²¹⁰Pb method. A reference chronology from Sellafield-derived radionuclides serves to assess the ²¹⁰Pb dating models. The study uses the mass depth scale and solves the estimation of the unsupported fraction of ²¹⁰Pb needed for the models. The profile was very irregular, and a cluster analysis led to an ambiguous use of the piecewise CFCS model. The inventory was incomplete and in an unsteady state, but the CRS model can be tentatively applied with the reference SAR and the reference date methods, although also with ambiguous results. TERESA can explicitly handle ²¹⁰Pb_exc fluxes and sedimentation rates that vary over time and shows the best performance with insightful outputs.

1. Introduction

1.1. The Fundamentals of the ²¹⁰Pb-Based Method for Dating Recent Sediments

The ²¹⁰Pb dating method is widely used to date recent sediments (<100–150 yr), as it can provide an absolute age determination. The method was first proposed for dating glacier ice [1]. Then, it was applied to lacustrine and marine sediments [2,3].

The technique is based on the geochemical cycling of ²¹⁰Pb in nature, and distinguishes two components in the total ²¹⁰Pb content of a sediment sample [4]: (i) the so-called supported ²¹⁰Pb, in secular equilibrium with its parent radionuclide ²²⁶Ra, and (ii) the unsupported or excess ²¹⁰Pb, labelled hereafter as ²¹⁰Pb_exc, from atmospheric meteoric origin. This applies to aquatic sediments that remain permanently saturated, so that the diffusion and scape of ²²²Rn is generally negligible.

Meteoric ²¹⁰Pb comes from the decay of ²²²Rn exhaled from the continental crust and dispersed in the atmosphere. The path of meteoric ²¹⁰Pb before reaching the studied sediment can be very complex. In aquatic environments, it can be bound to solid particles at different stages of this path, and it can spend some time in transient deposits before reaching the site studied [5].

The ²¹⁰Pb_exc is ultimately a meteoric ²¹⁰Pb, but for radiometric dating purposes, its history begins when it is incorporated into the sediment–water interface (SWI). The ²¹⁰Pb_exc content in the accreting sediment at the SWI is called the initial activity concentration, $A_0$ (with typical units of Bq kg⁻¹). From then on, it undergoes radioactive decay (T_1/2 = 22.23 ± 0.12 yr for ²¹⁰Pb), while the former layer in the SWI is buried by accreting materials. Post-depositional redistribution processes may be active or not, depending on the particular sedimentary scenario.

1.2. ²¹⁰Pb-Based Dating Models for Recent Sediments

A capital question is that the empirical depth profile of mass-activity concentrations of ²¹⁰Pb_exc in the sediment does not unambiguously define a chronology (an age-depth relationship). The chronology cannot be decoded without a series of assumptions on the prevailing sedimentary conditions (e.g., on the temporal variability in accretion rates and $A_0$, and on the existence and parameterization of mixing, diffusion, etc.). There are different possibilities for selecting such assumptions, each functional choice is referred to as a ²¹⁰Pb-based dating model, and different models potentially lead to different results.

In recent decades, more than 20 models have been postulated to cope with the wide diversity of sedimentary conditions found in nature. This includes the piecewise versions of some models. Some reviews can be found, among others, in Mabit et al. (2014) [6] and Arias-Ortiz et al. (2018) [7]. A summary of the set of models used in this work is presented in Supplementary S1.

The unifying theory considers the sediment as a continuous medium, where the mass conservation of a particle-bound tracer is given by the classical advection-diffusion equation, with the continuity of fluxes at the SWI as a boundary condition [8]. In the absence of post-depositional processes, this condition implies that new inputs of tracers and matter are ideally deposited over the previously existing material and mathematically expressed as follows.

(1)$F=A_{0}w$

with $F$ being the ²¹⁰Pb_exc flux (with typical units of Bq m⁻²yr⁻¹), and $w$ being the sediment accumulation rate (SAR) (with typical units of kg m⁻²yr⁻¹ or g cm⁻² yr⁻¹). In some cases, Equation (1) does not hold, and the fast depth distribution of fluxes through kinetic reactive transport in a poly-phasic medium must be considered [9,10].

Due to the intrinsic ambiguity in the ²¹⁰Pb-based chronologies, an external validation is mandatory [11]. Independent chronostratigraphic markers, such as some bomb fallout radionuclides (¹³⁷Cs, ²⁴¹Am, and ^239+240Pu) are widely used for this goal [12,13,14,15]. Data for external validation are not always available and a large number of published studies only report the ²¹⁰Pb-based chronology.

In the context described above, the ²¹⁰Pb dating technique has been in a continuous process of construction and revision. However, practical applications to environmental studies have focused on the use of a reduced set of models, mainly the constant rate of supply (CRS) model, and the constant flux with constant sedimentation (CFCS) model [4,16]. Both models have a simple mathematical formulation, which partially explains their popularity. They assume Equation (1) with constant flux. Additionally, CFCS assumes a constant SAR, and CRS assumes a steady-state total inventory of ²¹⁰Pb_exc.

1.3. Temporal Variability in ²¹⁰Pb_exc Fluxes

The paradigm of constant flux has been almost universally accepted, including models with mixing or diffusion. This is a mathematically convenient assumption since it allows for a simple analytical formulation of the model. However, there is wide empirical evidence that the most common sedimentary conditions involve fluxes and SARs that largely vary over time and achieve a positive statistical correlation among them [5].

The above poses the question of studying the model errors that arise from a partial or null achievement of the model assumptions. These studies have shown that CFCS and CRS models can produce reasonable proxies to the true chronology when the temporal variability in ²¹⁰Pb_exc fluxes is randomly distributed in the time line, but these models fail with continuous trends of increase or decrease in fluxes [17,18].

The only dating models that explicitly consider $F$ varying with time are: (i) the Constant Initial Concentration (CIC) model that assumes a constant $A_0$ [1]; (ii) the Sediment Isotope Tomography (SIT) model, which claims its ability to work with simultaneous and independent temporal variability in F and w [19]; and (iii) Time Estimates with Random Entries of Sediment and Activities (TERESA), which relies on the independent variability in $A_0$ and SARs [18,20].

The SIT model is not considered here since it has been shown to lack a solid physical basis [21]. The CIC model has a simple mathematical formulation and is a useful tool for preliminary evaluations and hypothesis testing. It has shown its use in real applications, but with complex ²¹⁰Pb_exc versus depth profiles, this model leads to age reversals. Therefore, it is not considered in this study.

In the case of steeped changes in the mean sedimentary conditions, the piecewise versions of the CFCS, CRS, and (multimodal) TERESA can be used, but this latter is the only model that can handle continuous trends of change in ²¹⁰Pb_exc fluxes [17,18].

For the scientific community working with the radiometric dating of recent sediments, comparative studies on the performance of the different models in scenarios with some degree of complexity can be insightful. This is the aim of the present study.

1.4. A Challenging Sediment History for Testing ²¹⁰Pb-Based Dating Models

Kershaw et al. (1990) [22] sampled sediment cores in 1988 in Maryport. The site had been dredged in the early 1950s and dredged again shortly after sampling. In this period of time, the sediment grew more than 2 m without disturbance, keeping a precious high-resolution record of a wide set of radioactive tracers, preserving the history of Sellafield discharges and the fingerprint of the wastes released by the phosphate plant at Whitehaven, with enhanced levels of ²²⁶Ra and ²¹⁰Pb. The authors were unable to apply the ²¹⁰Pb method, leaving an open problem that is faced here.

This sedimentary system is particularly challenging due to the following: (i) there is no total inventory of ²¹⁰Pb_exc due to dredging and incomplete recovery; (ii) the constant equivalent flux required to conform the measured inventory is about two orders of magnitude higher than the atmospheric deposition in the area (details are reported in this work); (iii) the ²¹⁰Pb_exc profile was highly irregular; (iv) the temporal resolution defined by the slicing is of a few months, much higher than the typical one of a few years; (v) the ²²⁶Ra baseline is not uniform with depth, and it reflects the impact of the wastes of the phosphate industry released to the Irish Sea in Whitehaven; (vi) ²¹⁰Pb is also present in such releases in radioactive disequilibrium with ²²⁶Ra; and (vii) ²¹⁰Pb_exc fluxes may not be constant because of natural variability and the impacts of the nearby phosphate industry.

The first point breaks one of the basic assumptions of the CRS model. In the 1990s, the use of ²¹⁰Pb dating methods was dominated by the view of ²¹⁰Pb_exc fluxes identified with constant atmospheric deposition. Thus, the criteria for the reliability of the CRS and CFCS assumptions included that the measured inventory should be within a range expected from regional atmospheric fallout rate and that multiple cores within a lake, but from areas of different sediment accumulation rates, should have similar ²¹⁰Pb_exc inventories. Therefore, the above point (ii) was particularly challenging. Point (iii) limits the use of the CFCS model. The remaining points concern all of the models that require constant flux.

Kershaw et al. (1990) [22] provided a reliable age-depth model based on the anthropogenic radionuclides released by Sellafield. Therefore, their data set gives a unique opportunity to test the performance of the ²¹⁰Pb dating models with lessons learnt in the last three decades. The evidence reported by the authors allows discarding the family of models considering mixing or diffusion, so this paper focusses on the use of the CFCS, CRS, and TERESA models, by applying up-to-date methods.

2. Materials and Methods

2.1. The Empirical Data Set

This work uses a compilation of literature data including (i) the reported physical and radiometric measurements in the sediment core sampled in Maryport [22]; (ii) the records of liquid radioactive discharges in Sellafield (data from Gray et al., 1995 [23]); and (iii) the releases of natural radionuclides from the phosphate plant at Whitehaven (data from Camplin et al., 1996 [24]). Figure S1, in the Supplementary Materials shows the locations referred to above along with the sampling site of the core studied. A brief summary of the data is presented below.

The sediment studied was sampled in February 1988 at Senhouse Dock, Maryport, at a site that had been dredged in the early 1950s. Two cores were sampled 1 m apart, labelled M1 and M2, with lengths of 210 cm and 108 cm, respectively. The reported shortening during the coring operation was about 20%. Kershaw et al. (1990) [22] reported mass-activity concentrations versus depth profiles for 21 radionuclides by merging sections of the two cores. They used the following sequences (with the range of depths given in brakes, with units of cm):

• CS1: M1(0–6) + M2(6–108) + M1(100–210), for dry/wet weight quotients, percentage of sands, LOI, ⁴⁰K, ¹⁵⁴Eu, ¹⁵⁵Eu, ¹²⁵Sb, ⁶⁰Co, ¹⁰⁶Ru, ¹⁴⁴Ce, ¹³⁴Cs, and ¹³⁷Cs.

• CS2: M1(0–6) + M2(6–60) + M1(60–210), for ²²⁶Ra, ²¹⁰Pb_exc, ²³⁸U, and ²³⁴Th.

• CS3: M2(0–108) + M1(100–210), for ^239,240Pu, ²³⁸Pu, ²⁴¹Am, ²⁴⁴Cm, and ²⁴²Cm.

The present work uses ⁴⁰K, ²¹⁰Pb, ²²⁶Ra, ²³⁸Pu, ^239,240Pu, ¹³⁷Cs, and ¹³⁴Cs. Raw data will be presented in tables and graphs throughout the article, along with results from the modelling study. The authors reported an age–depth model based on the interpretation of fingerprints of the artificial radionuclides in terms of the history of their releases at Sellafield. It will be used in this work for the sake of comparison.

Kershaw et al. (1990) [22] reported the ratio between dry weight and wet weight, labelled q, for each sediment section. Assuming a three-phase medium (solids, pore fluid, and air pores) of uniform density and volume additivity, and neglecting corrections by precipitated salts, the dry bulk density, ${\rho }_b$, can be estimated as

(2)${\rho }_b={\displaystyle \frac{1-{\varphi }_{air}}{\frac{1}{{\rho }_s}+\frac{1-q}{q{\rho }_w}}}$

where ${\varphi }_{air}$ is the fraction of the bulk volume occupied by air pores, ${\rho }_s$ is the density of solids (typically ~2.5 g cm⁻³), and ${\rho }_w$ is the density of the saline pore fluid (taken as 1.02 g cm⁻³). For a water-saturated sediment, ${\varphi }_{air}$ = 0, but higher values cannot be discarded in the upper sediment layers since the sampling was carried out during low water, with the sediment exposed to the atmosphere.

It was not clear whether the shortening (of ~20%) was due to material loss at the base of the core prior to closing the spring-loaded doors, or to the extrusion of part of the pore fluid, or both. Kershaw et al. (1990) [22] observed negligible bioturbation. The X-ray radiographs revealed a well-preserved sedimentary structure; thus, it can be ensured that the mass depth, $m$ [M L⁻²], remains invariant and can be used as the natural depth scale for dating models. This is consistent with the fact that bulk density (Equation (2)) represents the state of the sediment at the time of the slicing, which may differ from that of the natural conditions. The mass depth accumulated from the sediment–water interface (SWI) to a depth z (at the time of slicing), $m\left(z\right)$, is

(3)$m\left(z\right)={{\displaystyle \int }}_0^z{\rho }_b\left(z^{\prime }\right)dz\prime$

Table S1, in Supplementary Materials lists the q values as reported by Kershaw et al. (1990) [22], and computes the bulk densities by Equation (2). Figure 1 plots ${\rho }_b\left(z\right)$ along with the depth profiles of the normalised concentrations of ⁴⁰K and the sand fraction (particles with diameters > 63 µm). Normalisation here is defined as the ratio between the value of the slice and its arithmetic mean for all of the measured slices in the core. It is used only to facilitate comparison in the same plot of magnitudes with different units and scales [15]. The insights that Figure 1 brings are discussed in the Section 3.

Table S2, in the Supplementary Materials summarises the raw data for ²²⁶Ra and ²¹⁰Pb_exc reported by Kershaw et al. (1990) [22]. Figure 2 shows the ²²⁶Ra content in the sediment core, along with the historical annual liquid releases from the phosphate plant in Whitehaven, Cumbria, as estimated by Camplin et al. (1996) [24]. These releases, discharged to the sea via a pipeline, consisted of a phosphogypsum (PG) slurry and acid wastes from the further processing of the phosphoric acid. The data for the ²¹⁰Pb releases are not corrected by radioactive decay and ingrowth from ²²⁶Ra. Before 1987, these discharges were best estimate values and were not based on measurements at the pipeline.

2.2. ²¹⁰Pb-Based Radiometric Dating Models

In the sediment studied, along with the natural baseline values of ²²⁶Ra, there is a contribution of anthropogenic ²²⁶Ra (from the releases to the Irish Sea of the phosphate plant at Whitehaven). The latter is mostly associated with PG particles [24], and it is in a noticeable radioactive disequilibrium with its progeny radionuclide ²¹⁰Pb (²²⁶Ra activity is about a factor of 1.6 higher than that of ²¹⁰Pb at the time of the releases, as seen in Figure 2). The sketch in Figure S2, in the Supplementary Materials, illustrates the different components of the ²²⁶Ra and ²¹⁰Pb mass-activity concentrations in the material reaching the SWI. The deficit of supported ²¹⁰Pb (or the excess of ²²⁶Ra) is compensated for by a fraction of equal size from meteoric ²¹⁰Pb. The radioactive decay of this fraction of meteoric ²¹⁰Pb is compensated for by the ingrowth of ²¹⁰Pb from the excess ²²⁶Ra. Consequently, the total ²²⁶Ra content must be subtracted from the total ²¹⁰Pb to obtain the initial activity concentration $A_0$ to be used in dating models.

Along with the liquid waste released to the Irish Sea, the phosphate plant also delivers to the atmosphere ²²²Rn gas and phosphate dust generated during the loading and processing of phosphate ore [25]. The former can be dispersed over long distances, contributing to the enhancement of the atmospheric deposition of ²¹⁰Pb.

Kershaw et al. (1990) [22] discarded bioturbation, the core had a well-preserved sedimentary structure, and their age–depth model was based upon the linear relationship among the depth profiles of Sellafield-derived radionuclides in the core and the history of their releases. Consequently, only the family of ²¹⁰Pb-based models that excludes post-depositional redistribution is considered here. It comprises the models CFCS, CRS, and TERESA. For the CRS model, potential insights from the application of the reference date and the reference SAR methods are explored in this study [17,26]. For the models studied, a summary of their constitutive equations is presented as Supplementary S1, in the Supplementary Materials and full details can be found in the references provided above. For TERESA, the codes for numerical solutions are available as Supplementary Materials in Abril (2020) [18]. Bruel and Sabatier (2020) [27] presented an R-based software package, serac (https://github.com/rosalieb/serac), which allows the application of the CRS, CIC, and CFCS models, including their piecewise versions.

3. Results and Discussion

3.1. The Physical Porous Medium

The bulk density (Figure 1) was calculated assuming saturation conditions. It shows a strong departure from the typical profiles of natural compaction with an asymptotical trend of increasing with depth [8]. This feature can indicate (i) partial saturation conditions, with non-null values of ${\varphi }_{air}$ in the upper slices; (ii) core shortening with loss of interstitial water during coring and slicing; or (iii) disruptions of natural self-compaction due to changes in the sedimentary conditions.

The first factor cannot be discarded for the uppermost slices, since the sampling was carried out during low tide, with sediments exposed to the atmosphere. However, and due to the low infiltration velocities (of a few mm/hour; [28]), it can hardly correct the empirical ${\rho }_b\left(z\right)$ in the 0–30 cm to fit the asymptotical increase pattern. The lack of empirical data on ${\varphi }_{air}$ does not allow for reliable corrections. The present discussion serves to highlight the interest of determining ${\rho }_b$ from the dry weight of a known bulk volume. Then, ${\varphi }_{air}$ can be estimated by comparing the empirical ${\rho }_b$ against the one estimated by the three-phase model (Equation (2)). Unlike the other factors above, air pores affect the estimation of the bulk volume under the slicing conditions, and then, the estimation of the mass depth scale.

The second factor potentially affecting ${\rho }_b\left(z\right)$ seems reasonable, but of moderate extent due to the agreement in the geochemical analyses of interstitial water reported by Kershaw et al. (1990) [22] for this core and for a supplementary one, of 20 cm length, sampled at the same site. The third factor would be supported by the observed drastic change in the sand fraction in the upper 35 cm of the core (Figure 1), with the limitation of the available partial data.

The mass-activity concentrations of ⁴⁰K (Figure 1) are an indicator of the clay content of the sediment. They show relatively uniform values throughout the core, with a relative standard deviation of 8.0%.

3.2. Natural and Anthropogenic Sources of ²²⁶Ra and ²¹⁰Pb

The total of ²²⁶Ra is the relevant magnitude for estimating the unsupported fraction of ²¹⁰Pb useful for radiometric dating. The aim here is to provide an approximate estimation of the different components of ²²⁶Ra in this core, affected by the ²²⁶Ra releases from the phosphate plant at Whitehaven.

The ²²⁶Ra data in the core (Figure 2) have a poor spatial resolution below a 40 cm depth, but roughly follow a pattern of increase, in parallel with the historical releases. However, this covariation breaks down in the upper ~45 cm of the core, probably reflecting changes in local sedimentary conditions. This point will be discussed further below in light of the results of the dating models.

The ²²⁶Ra values in the sediment core can be seen as a natural background, enhanced by the waste from the phosphate plant that reaches the sampling site, most likely as PG particles. The background component can be roughly estimated from the average of the four deepest data points, being of ~15 Bq kg⁻¹. The maximum value observed for anthropogenic ²²⁶Ra (above the baseline values) is 21.1 Bq kg⁻¹ at $m=$ 48.5 g cm⁻². The total inventory of anthropogenic ²²⁶Ra can be estimated as ~16.5 × 10³ Bq m⁻², using linear interpolation for missing measurements. This figure can be compared with the integrated discharges from 1954 to 1987, both included, being of 16.6 × 10¹² Bq. Their ratio, $P_R$ ~1.0 × 10⁻⁹ m⁻², can be interpreted as the probability that a unit of ²²⁶Ra released at Whitehaven reached an area of 1.0 m² at the sampling site. Note that corrections by radioactive decay are negligible for ²²⁶Ra (T_1/2 = 1600 yr) for the involved time scales of a few decades.

If assumed constant over time, $P_R$ can be used to convert data on ²²⁶Ra releases, $D_{\mathrm{Ra}}\left(t\right)$, into fluxes on anthropogenic ²²⁶Ra onto the SWI: $F_{Ra}\left(t\right)=P_R{D}_{Ra}\left(t\right)$. Finally, the last can be converted into mass-activity concentrations in the sediment by dividing by the sedimentation rate, $w$ [M L⁻² T⁻¹]: $A_{Ra}=F_{Ra}/w$. To obtain the total ²²⁶Ra concentration, the natural background level must be added.

It can be assumed that ²¹⁰Pb released at Whitehaven accompanies ²²⁶Ra in the PG particles. Due to the strong disequilibrium between ²²⁶Ra and ²¹⁰Pb in the releases (Figure 2), ²²⁶Ra would reach the sediment core in excess with respect to companion ²¹⁰Pb. Note that this is a different concept from the above excess with respect to the ²²⁶Ra baseline values. The size of such a disequilibrium can be estimated as follows. The fraction of excess ²²⁶Ra at the date of the release, $t_r$, can be estimated from the data for discharges in Figure 2, and converted into fluxes into the sediment, as above, and finally into the mass-activity concentration of ²²⁶Ra that remains in excess with respect to the Whitehaven derived ²¹⁰Pb: $A_{Ra}^{exc}={\displaystyle \frac{P_R}{w}} \left[D_{Ra}\left(t_r\right)-D_{Pb}\left(t_r\right)\right]e^{-\lambda \left(t_s-t_r\right)}$. For a constant $w$, such $A_{Ra}^{exc}$ is found at mass depth $m=w·\left(t_s-t_r\right)$.

The results of the above estimates are shown in Figure 3. They have been calculated with a constant sedimentation rate of $w$ = 4.5 g cm⁻²yr⁻¹. The model does not account for the likely changes in local sedimentary conditions that affect the upper ~40 cm of the core, but provides a reasonable estimate of the relative size of the various components in the mass-activity concentrations of ²²⁶Ra. To the extent that ²²⁶Ra and ²¹⁰Pb released in Whitehaven follow a similar dispersion path in the sea (what is expected, at least for the PG slurry), such discharges are a source of total ²¹⁰Pb in the sediment core, but not of ²¹⁰Pb_exc.

It is worth comparing the above analysis with that reported for the same core by McCartney et al. (1990) [29]. These authors assumed a lower baseline value of ~10 Bq kg⁻¹ for ²²⁶Ra and compared the estimated anthropogenic contribution with the records of U releases from Whitehaven, which were the only data available at that time. They assumed secular equilibrium in the U decay chain, which is not the case, as later shown by Camplin et al. (1996) [24]. Comparison was done in terms of relative variations around their respective mean values and using the chronology reported by Kershaw et al. (1990) [22] (presented in the next subsection). The discrepancies in the upper core were less apparent than in Figure 3, but the main reason is the different behaviour of the ²²⁶Ra and ²³⁸U discharges. Thus, while ²³⁸U discharges show a noticeable relative minimum in 1979–1981, ²²⁶Ra and ²¹⁰Pb releases increased in the same period (Figure 2).

3.3. Sellafield-Derived Radionuclides and the Reference Age Model

Kershaw et al. (1990) [22] presented an age-depth model for the 0–140 cm of the core based on the comparison of radionuclide profiles in the core with their decay-corrected Sellafield discharges. The authors used time marks for ²⁴¹Am, ^239,240Pu, ¹⁰⁶Ru, ¹³⁷Cs, ¹³⁴Cs, ⁶⁰Co and the ratios ¹³⁷Cs/¹³⁴Cs, ¹³⁷Cs/Pu(α), ¹⁵⁴Eu/¹⁵⁵Eu, and ^239,240Pu/²³⁸Pu. They considered a lag time from discharges in Sellafield to sedimentation in Maryport, which was specific for each radionuclide. In this chronology, a period of one year corresponded to depth intervals ranging from 3 cm up to 10 cm, comprising an integer number of sediment slices. The chronology was reported without associated uncertainties, and it is summarised in Table S3, in the Supplementary Materials.

The chronology has been translated to the mass depth scale (Equation (3)), which allows estimating SAR records. They are shown in Figure 4, along with the comparison of the core profiles of ^239,240Pu and ²³⁸Pu with those estimated from the history of releases. The latter involve the age model, a lag time of 0.5 yr, and a fitted scaling factor of 71 × 10⁻¹² yr · kg⁻¹, which is the mean ratio $P_R/w$, with $P_R$ in the range (3.5–4.0) × 10⁻⁹ m⁻² for these isotopes of Pu. Note that the data in the core have a resolution higher than the data from annual discharges. Despite these limitations and the oversimplified estimation of fluxes onto the SWI, there is general reasonable agreement, particularly in the position of the major peaks. This supports the assumption of negligible post-depositional redistribution of highly particle-reactive radionuclides (such as Pu and Pb isotopes), as well as the use of the age–mass depth model as a reasonable proxy to the true chronology.

The comparison between the measured profile and the decay-corrected discharge records for ¹³⁷Cs and ¹³⁴Cs is shown in Figure S3, in the Supplementary Materials. The age–mass depth model is used with a lag time of 1.5 years and a global scaling factor of 2.6 × 10⁻¹² yr · kg⁻¹. There is an overall acceptable agreement, but the description of ¹³⁴Cs in the upper region of the core is particularly poor. Discrepancies are also apparent for ¹³⁷Cs in the deeper part of its profile.

Several limitations can be posed to the chronology by Kershaw et al. (1990) [22]: (i) The identification in the measured profiles of singular fingerprints in the history of the releases is to some extent subjective, as noticed by the authors, and they can be masked by natural variability in the dynamics of marine dispersion and the sedimentary conditions. (ii) The authors used different sequences combining samples from M1 and M2 cores for different radiotracers (see the Section 2) that may have been differently affected by shortening during coring. (iii) The need to use different time lags for each radiotracer may mask pre- or post-depositional processes, and they affect the estimation of ages and SARs. (iv) The use of an integer number of sediment layers to define a period of one year may introduce some bias. (v) SARs must be considered as semiquantitative estimates, as they can be largely affected by propagated uncertainties (not reported in the age model) and model errors [30].

Despite the limitations above, the combined use of a wide set of radiotracers and isotopic ratios strengthens the global chronology, which can be considered as a reliable proxy to the true one. It will serve as a reference for comparison with the ²¹⁰Pb-based chronologies.

3.4. ²¹⁰Pb-Based Dating Models

3.4.1. Inventory of ²¹⁰Pb_exc

Data on ²¹⁰Pb_exc reported by Kershaw et al. (1990) [22] are summarised in Table S2 (Supplementary Materials), and plotted in Figure 5. The inventory of ²¹⁰Pb_exc in the recovered core, estimated with the bulk density from Figure 1, is ${\mathsf{\Sigma }}_C=$ 118.9 ± 1.1 kBq m⁻². In the deepest sections of the core, ²¹⁰Pb_exc is well above 40 Bq kg⁻¹. Therefore, it can be inferred that the core did not reach the level of dredging in the early 1950s, so the above inventory was accumulated in a time lapse $t_C$ < 38 years (the precise date of the dredging was not reported). This allows estimating a lower limit for the equivalent constant flux at the SWI: $F_{eq}>\lambda {\mathsf{\Sigma }}_C/\left(1-e^{-\lambda t_C}\right)$. This yields $F_{eq}>$ 5330 ± 50 Bq m⁻²yr⁻¹. This value is about two orders of magnitude higher than the atmospheric flux of 85 Bq m⁻²yr⁻¹ reported for this region [22]. Moreover, McCartney et al. (1990) [29] reported ²¹⁰Pb_exc inventories in 28 sediment cores sampled in the eastern Irish Sea in the period 1982–1987. They ranged from 1.0 up to 14.1 kBq m⁻², being one or two orders of magnitude lower than the value found in the studied core. This was a challenging context for the application of ²¹⁰Pb dating models, since at those dates the common view identified ²¹⁰Pb_exc fluxes with constant atmospheric deposition [29,31]. Therefore, Kershaw et al. (1990) [22] concluded that the flux of ²¹⁰Pb_exc appeared to be enhanced by a local anthropogenic source and this posed a limitation on the use of the ²¹⁰Pb_exc dating technique in the studied area.

It has been shown (Figure 3) that the releases to the sea from the phosphate plant at Whitehaven are enhancing the levels of ²²⁶Ra and total ²¹⁰Pb in the Maryport core, but they are not necessarily a source of ²¹⁰Pb_exc for this sediment. The idea of a potential enhancement of ²¹⁰Pb_exc fluxes through preferential (with respect to ²²⁶Ra) maritime transport and to atmospheric releases will be discussed further, in the light of the results of the TERESA model. It is worth noting that the values of ²¹⁰Pb_exc in the upper region of the core are <110 Bq kg⁻¹, which are common values for sediments without anthropogenic input. In particular, the figure above compares well with the values reported for saltmarsh sediments in the Lynher-Tamar estuary, UK [32].

Abril and Brunskill (2014) [5] reported empirical evidence that a positive correlation between ²¹⁰Pb_exc fluxes and SARs was a widespread sedimentary condition. In the case studied, the harbour acts as an efficient sediment trap, receiving new inputs of suspended loads with each tidal cycle and capturing a large fraction of them by settling. The retrieved core has a mass depth of 157.6 g cm⁻² that has been accumulated in $t_C$ < 38 years, which yields a mean SAR > 4.15 g cm⁻²yr⁻¹. The matter particles that reach the SWI carry meteoric ²¹⁰Pb, bound to their surfaces at different stages of their path. Thus, the extraordinarily high value of $F_{eq}$ is consistent with the also extraordinarily high value of SAR. Note that factors that affect the amount and provenance of the mass flows (such as dredging and water works in the area) can produce changes in the ²¹⁰Pb_exc fluxes.

3.4.2. The Constant Flux with Constant Sedimentation Model

The CFCS model has been applied to the ²¹⁰Pb_exc versus mass depth profile, hereafter denoted as $A\left(m\right)$, including its piecewise version, which splits $A\left(m\right)$ into two and three transects. This is shown in Figure 5 that plots $\mathrm{Ln}\left[A\left(m\right)\right]$. The fitted parameters are reported in Table 1, along with the inferred sedimentary conditions (the model only captures the mean values of $F, A_0$, and $w$). They are reported in Table 1. The so derived chronologies are depicted in the second panel of Figure 5, along with the reference chronology, adapted from Kershaw et al. (1990) [22].

The dispersion of data in transect 1 opens the question of splitting the transect into several parts. In transects 2 and 3, the data density is lower and the separation of both transects is mostly supported by the influential point at m = 84.5 g cm⁻². The alternative interpretation of a joined transect 2 + 3 results in a better agreement with the reference chronology (Figure 5). In this two-transect version, the two regions are separated by a jump discontinuity, reflecting steeped changes in $A_0$. Possible changes in SARs cannot be resolved from the uncertainties involved (Table 1). A more detailed discussion will be presented below in light of the results from the TERESA model.

All fits show moderately strong correlations, with p values below 0.002 (Table 1), but this is not a proof of the performance of the model, as shown by the failure of the chronologies, particularly from the global CFCS and the piecewise version with three transects. Such a failure indicates that the model assumption of a constant mean value of fluxes and SARs in the transect (over which temporal variability can be randomly distributed) does not hold for the entire core (1-transect CFCS) nor for the three transects ${\mathrm{r}}_1, {\mathrm{r}}_2$, and ${\mathrm{r}}_3$. This is referred to as a model error. However, the two-transect approach seems to capture the gross features of the prevailing sedimentary conditions.

3.4.3. The Constant Rate of Supply Model

In the Maryport core, the total inventory (see the mathematical definition in Supplementary S1) did not reach a steady state because of dredging in the early 1950s. Additionally, the retrieval was incomplete. However, as the CRS model assumes non-post-depositional redistribution, the previous history does not affect the development of the $A\left(m\right)$ profile. Thus, the extrapolation of the $A\left(m\right)$ profile using the SAR value given by the trend of exponential decrease in the deepest measured slices allows for a gross estimate of the missing inventory that would have existed under steady-state conditions and constant flux [12]. This allows for preliminary dating. That is, the reference SAR method does not distinguish if the missing inventory actually exists or not.

When a reference date is available, the fraction of the inventory lying from the SWI to the depth of the known date can be used for estimating the equivalent constant flux. Under the assumption that such a flux has remained constant over time (or with randomly distributed temporal variability), it can be used to estimate the equivalent steady-state total inventory. This method requires external data (the reference date). The so-generated CRS chronology necessarily matches the reference date, and it can be considered reliable for younger ages as long as the assumption of constant flux holds. For older ages, the method applies an extrapolation that must be handled with care.

The CRS model has been applied to the sediment core from Maryport. Interpolations were used for missing data within the measured slices. The measured part of the total inventory is 118.9 ± 1.1 kBq m⁻². The reference SAR method was used with the reference SAR for the global profile, the combined transect 2 + 3, and the third transect (r_G; r₂₊₃ and r₃ in Figure 5, respectively), with the fitted parameters reported in Table 1. This produces a ‘missed’ inventory of 122 ± 11, 61 ± 6, and 41 ± 6 kBq m⁻² for r_G, r₂₊₃, and r₃, respectively. The resulting chronologies are shown in Figure 6. The one based on r₂₊₃ gives a reasonably good proxy to the reference chronology, while the others are divergent. In all cases, the ages are underestimated for the younger slices.

For applying the reference date method, the age of the peak of ²³⁸Pu at m = 67.4 g cm⁻² (see Figure 4) is known to be 15 yr (including the 0.5 yr lag time). A second application of the method uses the first pronounced peak at m = 40.1 g cm⁻² with a known age of 9 yr. Other time marks could have been chosen, but the present selection serves to illustrate the performance of this method. These reference dates allow us to estimate equivalent constant fluxes of 5250 ± 60 and 4640 ± 50 Bq m⁻²yr⁻¹ for reference dates of 15 and 9 years, respectively. Such fluxes lead to steady-state total inventories of 168.8 ± 2.0 and 149.4 ± 1.8 kBq m⁻², respectively. The resulting chronologies are shown in Figure 6. They intercept the reference chronology at the stated reference dates, but they largely deviate for deeper slices. The estimate based on the reference age of 15 years shows an overall acceptable agreement. Note that this method requires external input data (the reference date).

The reference SAR and reference date methods are relatively simple to apply, allowing for the construction of CRS chronologies that are sensitive to the reference values used. The above results will be better understood in light of the TERESA model.

3.4.4. TERESA Model

The TERESA model does not require the total inventories of ²¹⁰Pb_exc, but it needs continuous measurements. Interpolations have been used to replace missing data within the measured slices. Application cases can be found, among others, in [18,20,32,33,34].

Figure 7 compares the empirical ²¹⁰Pb_exc profile with the best fit produced by the TERESA model (the fitted parameters are reported in Table 2). Despite the large irregularities in the empirical profile, the model is overall in very good agreement with the data, but for the outlier at m = 84.5 g cm⁻². The model chronology is shown in Figure 7. It compares well with the reference ages, with a maximum deviation of about 1.5 yr.

The TERESA model explored about 10⁵ combinations of normal distributions for $A_0$ and SARs, each one containing 63 pairs of data whose best sorting downcore is numerically solved. Note that although the 63 values of $A_0$ and SARs are normally distributed around their mean values, their best sorting in the time line (downcore) can result in a random distribution or in a well-defined pattern of change over time. The use of normal distributions is a simplifying assumption that can be overcome with the multimodal version of the model. From all of the above possible combinations, in the studied core, the best performance is achieved with a relatively uniform value of SAR, and the palaeo-records of $A_0$ and fluxes are reported in Figure 8. This corresponds to the minimum of the $\chi$-function (Equation (S5) in the Supplementary Materials), as it is shown in Figure S4 (Supplementary Materials). The model parameters ${\overline{A}}_0$, $\overline{w}$, and s_A are well constrained, while the parametric lines for s_w define on the $\chi$-hypersurface an elongated valley in the range of 0.02–0.18, although with a relative minimum around 0.05 (Figure S4). The error estimate in the fitting parameters is based on the curvature of the parametric lines, following Bevington and Robinson (2003) [35]. The propagated uncertainties in the chronology are small in this case (<6.5%, average 0.9%, at 95% confidence level). They have been computed following the methods in Abril (2016) [20]. It is worth noting that, in addition to propagated uncertainties, in all of the models there are always nonquantifiable model errors [10]; in particular, the TERESA model uses normal distributions for $A_0$ and SARs that are only a proxy to the true (but unknown) ones.

3.4.5. Inter-Comparison of Models and Interpretation of Data

The global picture provided in Figure 8 is consistent with the cluster analysis of Table 1. The approach of two transects ($r_1$ and $r_{2+3}$), where discontinuity corresponds to a steep change in the mean value of $A_0$ (around which temporal variability is randomly distributed), is in excellent agreement with the two transects depicted in Figure 8. It is worth noting that in Table 1, the estimated mean values of SAR for the two transects are not different with statistical significance when accounting for the associated fitting errors.

Splitting the $r_{2+3}$ transect into two regions separated by the influential point at m = 84.5 g cm⁻², results in fluxes showing a slight trend of increase in both transects, which explains the failure of the piecewise CFCS model (i.e., the increasing $F$ is misinterpreted by the CFCS model as a lower $w$ value; see [18]).

When considering a single transect, the higher density of data in region $r_1$ biases the fit, so that SAR is overestimated.

From Figure 8, the mean value of $F$ weighted by the mass thickness is 5690 Bq m²yr⁻¹. For the application of the CRS model, the estimates by the reference SAR with $r_{2+3}$ and the reference age of 15 yr, led to constant fluxes of 5590 and 5250 Bq m²yr⁻¹, respectively, which are close to the above figure. However, this figure is an overestimation of the true fluxes in region $r_1$ and an underestimation in region $r_{2+3}$, which translates into CRS ages that are underestimated and overestimated in such regions, so that they intercept the reference chronology at intermediate points.

Figure 8 allows us to interpret the changes over time in sedimentary conditions. The mass flow that reaches the SWI at any time is a composite of materials with different provenances. That of marine origin may carry on an important fraction of the fluxes of the particle-reactive isotopes released by Sellafield, as well as enhanced levels of ²²⁶Ra and ²¹⁰Pb due to the releases from the phosphate plant at Whitehaven. This material also contains meteoric ²¹⁰Pb. The fluxes of terrigenous materials are initially free of anthropogenic radionuclides, although they can uptake them from the dissolved phase before settling at the SWI. Terrigenous materials have their particular signature of ²²⁶Ra and meteoric ²¹⁰Pb. The relative proportion of material with different provenances may change over time due to varying riverine and coastal inputs, urban development, dredging, and other waterworks conducted in the area.

A major change would have occurred around a mass depth of 37.5 g cm⁻² ($z ~$ 45 cm), dated around 1979 (Figure 8). Although the mean SAR would not have changed significantly, the relative contribution of terrigenous flows would have increased, displacing part of the mass flows with marine provenance. Supporting evidence includes the following: (i) ²²⁶Ra shows a pattern of increase with time (i.e., toward shallower depths in the sediment) in parallel to the history of releases in the phosphate plant (Figure 2), with the marine pathway being the only one possible for such a correlation. This pattern breaks at such a date, with ²²⁶Ra concentrations in sediments decreasing with time, while releases continue their trend of increase (Figure 2 and Figure 3). (ii) There is an increase in the proportion of sands in this upper region of the core (Figure 1), although the data are only partial and do not provide information on the mean grain sizes. (iii) There is a noticeable discontinuity in the ²¹⁰Pb_exc profile (Figure 7), which the TERESA model interprets as a stepped reduction in the mean value of $A_0$ to 75% of its previous value (Figure 8). The above characteristics are consistent with the fact that the ²¹⁰Pb_exc content in the terrigenous material is lower than in that of marine provenance. The exact reason for this stepped change in $A_0$ remains unidentified, but is likely linked to local changes in the Maryport area (e.g., the dredging of neighbouring areas, other waterworks, changes in the urban planning, etc.).

Below a 37.5 g cm⁻² mass depth, the major fingerprint in $A_0$ is the relative maximum at 84.5 g cm⁻² (dated around 1969; see Figure 7 and Figure 8), which co-occurs with a relative maximum in ²²⁶Ra. The likely reason would have been an exceptionally higher contribution of marine flows of matter. However, the global picture of fluxes and $A_0$ in this region can be interpreted as natural random variability around their respective constant mean values (Figure 8), as assumed in the $r_{2+3}$ approach in Figure 5. Adopting the approach of splitting this region into two transects, that of $r_3$ (see Figure 5) corresponds to a period of slight but continuous increase in both ²¹⁰Pb_exc fluxes and initial activity concentrations, which are statistically significant at the >99% confidence level. This can be interpreted as dynamic sedimentary conditions in which marine inputs progressively gained relevance. This trend is interrupted around 1969, reverting the contribution of marine inputs to previous stages, after which the transect $r_2$ re-encounters the trend of continuous and statistically significant increase in $A_0$. The reason for this stepped change in sedimentary conditions around 1969 remains unidentified, as does that of 1979. Both produced similar qualitative effects, although with less intensity in 1969.

From the analysis in Figure 3, the amount of meteoric ²¹⁰Pb needed to compensate for excess ²²⁶Ra in the releases of the phosphate plant is <6.0 Bq kg⁻¹, a negligible amount compared with the $A_0$ values depicted in Figure 8 (a range of 70–85 Bq kg⁻¹). However, the above $A_0$ values compensate for this effect. As commented above, the temporal variability in $A_0$ can be explained by changes in the relative proportion of marine and terrigeneous inputs. It is possible that periods of increasing $A_0$ could be partially contributed to by an enhancement of the ²¹⁰Pb_exc fluxes because of the gaseous emissions at Whitehaven. Similarly, preferential maritime transport (with respect to ²²⁶Ra) of the dissolved fraction of ²¹⁰Pb released to sea by the phosphate plant, could enhance the $A_0$ values at the core site. Preferential transport must be understood in terms of the interaction of these radionuclides with suspended matter and surficial sediments. However, although they cannot be discarded, if these contributions were the main reason for the patterns of increase in $A_0$, they could not explain the disturbances in 1969 and 1979.

It is worth noting the reasonably good performance of the TERESA model in this core, with a temporal resolution in the slicing of a few months. The fundamentals of the model rely on the independent variability in $A_0$ and SAR, supported by the wide empirical evidence provided in the work of Abril and Brunskill (2014) [5]. The last was from varved sediments, with typical temporal resolutions in the slicing of the order of one to a few years. In the studied core, the slicing is decoupled from any potential seasonal variability, so the variability in $A_0$ and SAR around their mean values remains essentially random and independent.

In Figure 5, Figure 6 and Figure 7, the reference chronology by Kershaw et al. (1990) [22], although only a reasonable proxy to the true (but unknown) one, can be used for independent validation of the chronologies based on ²¹⁰Pb. Without such an external validation, it is not possible to elucidate which of the different choices for transects for the piecewise CFCS is of the best performance (it has been shown that the p-value is not an objective criterion for model reliability). The application of the CRS model in this case is based on extrapolations for the estimation of the equivalent constant flux, with different possibilities for their choice that lead to different chronologies, so external validation is mandatory. The TERESA model involves more complex mathematics, but this is not a guarantee of its performance.

Independent validation must be, in all of the cases, an essential part in the methodology. However, this is not always possible, so methods to support the self-consistency of the ²¹⁰Pb-based chronologies are of interest. The above comparison of the TERESA model with the cluster analysis in Table 1 and the CRS projections (as shown in Figure S5, Supplementary Materials) provides a promising approach for this goal. This requires a good understanding of the involved model assumptions and how models perform when such assumptions are stressed (e.g., for the different types of temporal variability considered for the ²¹⁰Pb_exc fluxes). A holistic analysis of other magnitudes, such as bulk density, granulometry, or depth profiles of naturally occurring radionuclides, can support the ²¹⁰Pb-based chronologies.

4. Conclusions

The performance of the ²¹⁰Pb-based dating models CFCS, CRS, and TERESA has been tested in the challenging sedimentary scenario captured by a core sampled in 1988 in the Maryport harbour.

The Whitehaven phosphate plant delivers ²²⁶Ra and ²¹⁰Pb to the sea in radioactive disequilibrium, but the relevant magnitude for radiometric dating is the ²¹⁰Pb found in excess with respect to ²²⁶Ra at the time of its incorporation into the SWI. The analysis of the source term and the depth profiles of these radionuclides does not allow us to discard the applicability of ²¹⁰Pb-based models.

A cluster analysis of the $\mathrm{Ln}\left[A\left(m\right)\right]$ plot revealed potential stepped changes in fluxes, with jump and slope discontinuities. It supports the application of the piecewise version of the CFCS model, but with different possibilities for defining the set of clusters, which leads to different chronologies.

A total steady-state inventory is not defined for this core, but the CRS model can be applied in an exploratory way on the basis of extrapolations by using the well-known reference-SAR and reference-date methods. Different choices of the reference values lead to very different chronologies.

The TERESA model explicitly handles independent temporal variability in $A_0$ and SARs. It has shown the best performance. The model outputs historic records of $A_0$ and fluxes that allow for a better understanding of the limited performance of the CFCS and CRS models when their basic assumptions are stressed. It also serves to decode the gross features of changes in sedimentary conditions at the sampling site.

The temporal changes in $A_0$ and fluxes can be reasonably understood from the varying contributions of mass flows with marine and terrigeneous provenances, the former carrying on the ²²⁶Ra from the phosphate plant, and higher contents of excess ²¹⁰Pb.

The analysis of a sediment core by using the set of models, TERESA, CFCS, and CRS, along with the holistic study of other physic-chemical proxies, can strengthen the ²¹⁰Pb-based chronologies, when the potential and limitations of these models are properly understood.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are available in Supplementary Materials.

Conflicts of Interest

The author declares no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Bulk density in the Maryport sediment core, estimated by Equation (2) from the quotients of dry/wet weight reported by Kershaw et al. (1990) [22], and assuming saturation conditions. Also shown are the normalized values of 40K and the sand fraction (particles with diameter [greater than] 63 µm), with arithmetic means of 629.7 Bq kg−1 and 4.61%, respectively (raw data from the reference cited).

Figure 2. 226Ra in the Maryport sediment core (data from Kershaw et al., 1990 [22]). The error bars are smaller than the size of the markers. In the secondary axes, the annual liquid discharges of 226Ra and 210Pb from the phosphate plant at Whitehaven are presented (data from Camplin et al., 1996 [24]). The selected scales allow for comparison, but do not represent a mass-depth age model.

Figure 3. 226Ra in the Maryport sediment core, as in Figure 2. It is compared with the estimated concentrations using the model described in the text. It considers a natural baseline value (~15 Bq kg−1), and the anthropogenic contributions of the releases at Whitehaven, distinguishing for the latter the fraction that remains in excess with respect to 210Pb of the same origin. The model uses the discharge records from Figure 2, [Forumla omitted. See PDF.] = 1.0 × 10−9 m−2 and a constant [Forumla omitted. See PDF.] = 4.5 g cm−2yr−1.

Figure 4. 239,240Pu and 238Pu in the Maryport sediment core, compared to the records of their releases in Sellafield (data from Gray et al., 1995 [23]) by using the age model (plotted on the secondary Y-axis), a delay time of 0.5 yr and a global scaling factor of 71 × 10−12 yr · kg−1. Data for the core and chronology were adapted from Kershaw et al. (1990) [22]. The data for 238Pu are corrected by radioactive decay to the date of sampling. Vertical error bars are 1σ analytical uncertainties, and horizontal bars delimit the mass depth interval for the estimate of SARs.

Figure 4. 239,240Pu and 238Pu in the Maryport sediment core, compared to the records of their releases in Sellafield (data from Gray et al., 1995 [23]) by using the age model (plotted on the secondary Y-axis), a delay time of 0.5 yr and a global scaling factor of 71 × 10−12 yr · kg−1. Data for the core and chronology were adapted from Kershaw et al. (1990) [22]. The data for 238Pu are corrected by radioactive decay to the date of sampling. Vertical error bars are 1σ analytical uncertainties, and horizontal bars delimit the mass depth interval for the estimate of SARs.

Figure 5. Plot of [Forumla omitted. See PDF.] for the Maryport sediment core (upper panel), along with linear regressions applied to the global profile (rG, blue dotted line), to the three transects r1, r2, and r3, and to the aggregation of the two last (r2+3), identified with different colours (the red circle in the outline part of the gray and orange circles refers to r2+3). The fitted parameters are listed in Table 1. The (lower panel) depicts the resulting chronologies (with 1σ propagated error-bars) from the CFCS model and its piecewise version, along with the reference chronology (adapted from Kershaw et al., 1990 [22]). The piecewise CFCS model is applied with two (r1 plus r2+3) and three (r1, r2, and r3) transects.

Figure 6. Chronologies produced by the CRS model for the Maryport sediment core using the reference SAR and reference date methods. They are compared with the reference chronology of Kershaw et al. (1990) [22]. Propagated uncertainties are depicted only in two cases, for the sake of clarity.

Figure 7. Empirical 210Pbexc profile for the Maryport sediment core, and the best fit ([Forumla omitted. See PDF.] = 0.80) produced by the TERESA model using the parameters reported in Table 2. The model chronology is plotted in the secondary Y-axis, along with the reference by Kershaw et al. (1990) [22].

Figure 8. Records of [Forumla omitted. See PDF.] versus mass depth and fluxes versus model date (in secondary axis) produced by the TERESA model for the Maryport sediment core. Horizontal dashed lines are the mean values of fluxes for the transects above and below the 37.5 g cm−2 mass depth (of 4577 and 6079 Bq m−2yr−1, respectively).

Table 1

Cluster analysis of the Maryport sediment core.

Parameter (From Fit)		Sedimentary Conditions
A single transect (n = 44)
−s	0.0042 ± 0.0005	A 0	99.0 ± 3.6
b	4.5908 ± 0.0357	w	7.39 ± 0.93
Pearson r; p	−0.7726; 0.0000	F	7300 ± 950
Transect 1 (n = 24)
−s1	0.0079 ± 0.0022	A 0,1	99.6 ± 4.4
b 1	4.6001 ± 0.0441	w 1	3.94 ± 1.12
Pearson r; p	−0.6009; 0.0019	F 1	3900 ± 1100
m1 discontinuity	37.51
Transect 2 (n = 9)
−s2	0.0122 ± 0.0018	A 0,2	150 ± 14
b 2	4.7160 ± 0.0040	w 2	2.54 ± 0.38
Pearson r; p	−0.9299; 0.0003	F 2	3800 ± 670
m2 discontinuity	82.84
Transect 3 (n = 11)
−s3	0.0102 ± 0.0017	A 0,3	205 ± 29
b 3	4.4723 ± 0.0735	w 3	3.05 ± 0.75
Pearson r; p	−0.8923; 0.0002	F 3	6300 ± 1400
Transect 2 + 3 (n = 20)
−s2+3	0.0073 ± 0.0008	A 0,2+3	143 ± 14
B 2+3	4.6690 ± 0.0499	w 2+3	4.23 ± 0.44
Pearson r; p	−0.9132; 0.0000	F 2+3	6060 ± 870

In the Ln[A(m)] plot (Figure 5), A(m) is given in Bq kg⁻¹ and m in g cm⁻². A₀,_i and w_i (i = 1, 2, 3—it refers to transects or clusters) have units of Bq kg⁻¹ and g cm⁻²y⁻¹, respectively, and F_i are given in Bq m⁻²y⁻¹; s_i are slopes, and b_i the independent terms in the linear fit. The values and errors in A₀,_i for i > 1 have been estimated from Ln[A(m′)] plots with a translation to the origin (m′ = m − m₁ for i = 2; m′ = m − m₂ for i = 3). Full details on the methodology can be found in Abril (2019) [17].

Table 2

Entry parameters for the TERESA model and fitting-error estimates.

${\overline{\mathit{A}}}_0$ (Bq kg−1)	$\overline{\mathit{w}}$ (g cm−2 yr−1)	sA	sw
120.70 ± 0.16	4.540 ± 0.002	0.176 ± 0.001	0.050 ± 0.002

$\overline{A_0}$ and $\overline{w}$ are the arithmetic mean values of the initial concentration and SAR, respectively, and s_A, s_w are the corresponding relative standard deviations. The fit yields a χ value of 0.800.

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/jmse13010144/s1: Supplementary S1: A summary of the radiometric dating models; Figure S1: Map of the site; Figure S2: Sketch with the components of total ²¹⁰Pb; Figure S3: ¹³⁷Cs and ¹³⁴Cs in the sediment core; Figure S4: Mapping of the χ function for the TERESA model; Figure S5: Chronologies for the Maryport sediment core; Table S1: Bulk density; Table S2: ²¹⁰Pb data; Table S3: Chronology from Sellafield-derived radionuclides.

References

1. Goldberg, E.D. Geochronology with Pb-210. Proceedings of a Symposium of Radioactive Dating; International Atomic Energy Agency: Vienna, Austria, 1963; pp. 121-131.

2. Krishnaswamy, S.; Lal, D.; Martin, J.M.; Meybek, M. Geochronology of lake sediments. Earth Planet. Sci. Lett.; 1971; 11, pp. 407-414. [DOI: https://dx.doi.org/10.1016/0012-821X(71)90202-0]

3. Koide, M.; Soutar, A.; Goldberg, E.D. Marine geochronology with ²¹⁰Pb. Earth Planet. Sci. Lett.; 1972; 14, pp. 442-446.

4. Robbins, J.A. Geochemical and Geophysical applications of radioactive lead isotopes. Biochemistry of Lead in the Environment; Nriago, J.P. Elsevier: Amsterdam, The Netherlands, 1978; pp. 285-393.

5. Abril, J.M.; Brunskill, G.J. Evidence that excess ²¹⁰Pb flux varies with sediment accumulation rate and implications for dating recent sediments. J. Paleolimnol.; 2014; 52, pp. 121-137. [DOI: https://dx.doi.org/10.1007/s10933-014-9782-6]

6. Mabit, L.; Benmansour, M.; Abril, J.M.; Walling, D.E.; Meusburger, K.; Iurian, A.R.; Bernard, C.; Tarjan, S.; Owens, P.N.; Blake, W.H. et al. Fallout²¹⁰Pb as a soil and sediment tracer in catchment sediment budget investigations: A review. Earth Sci. Rev.; 2014; 138, pp. 335-351. [DOI: https://dx.doi.org/10.1016/j.earscirev.2014.06.007]

7. Arias-Ortiz, A.; Masqué, P.; Garcia-Orellana, J.; Serrano, O.; Mazarrasa, I.; Marbà, N.; Lovelock, C.E.; Lavery, P.S.; Duarte, C.M. Reviews and syntheses: ²¹⁰Pb-derived sediment and carbon accumulation rates in vegetated coastal ecosystems: Setting the record straight. Biogeoscience Discuss.; 2018; 15, pp. 6791-6818. [DOI: https://dx.doi.org/10.5194/bg-15-6791-2018]

8. Abril, J.M. A new theoretical treatment of compaction and the advective-diffusive processes in sediments. A reviewed basis for radiometric dating models. J. Paleolimnol.; 2003; 30, pp. 363-370. [DOI: https://dx.doi.org/10.1023/B:JOPL.0000007220.16908.d4]

9. Abril, J.M.; Gharbi, F. Radiometric dating of recent sediments: Beyond the boundary conditions. J. Paleolimnol.; 2012; 48, pp. 449-460. [DOI: https://dx.doi.org/10.1007/s10933-012-9622-5]

10. Abril, J.M.; Barros, H. Modelling the kinetic reactive transport of pollutants at the sediment-water interface. Applications with atmospheric fallout radionuclides. J. Environ. Radioact.; 2022; 242, 106790. [DOI: https://dx.doi.org/10.1016/j.jenvrad.2021.106790] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34890932]

11. Smith, J.N. Why should we believe ²¹⁰Pb sediment geochronologies?. J. Environ. Radioact.; 2001; 55, pp. 121-123. [DOI: https://dx.doi.org/10.1016/S0265-931X(00)00152-1]

12. Appleby, P.G. Dating recent sediments by 210Pb: Problems and solutions. Dating of Sediments and Determination of Sedimentation Rate; Illus, E. STUK A-145 STUK: Helsinki, Finland, 1998; pp. 7-24. Available online: https://www.osti.gov/etdeweb/servlets/purl/636061 (accessed on 11 January 2025).

13. Arnaud, F.; Magand, O.; Chapron, E.; Bertrand, S.; Boës, X.; Charlet, F.; Mélières, M.A. Radionuclide dating (²¹⁰Pb, ¹³⁷Cs, ²⁴¹Am) of recent lake sediments in a highly active geodynamic setting (Lakes Puyehue and Icalma—Chilean Lake District). Sci. Total Environ.; 2006; 366, pp. 837-850. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2005.08.013]

14. Tylmann, W.; Alicja Bonk, A.; Goslar, T.; Wulf, S.; Grosjean, M. Calibrating ²¹⁰Pb dating results with varve chronology and independent chronostratigraphic markers: Problems and implications. Quat. Geochronol.; 2016; 32, pp. 1-10. [DOI: https://dx.doi.org/10.1016/j.quageo.2015.11.004]

15. Abril, J.M.; San Miguel, E.G.; Ruiz-Cánovas, C.; Casas-Ruiz, M.; Bolívar, J.P. From floodplain to aquatic sediments: Radiogeochronological fingerprints in a sediment core from the mining impacted Sancho Reservoir (SW Spain). Sci. Total Environ.; 2018; 631, pp. 866-878. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.03.114] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29727997]

16. Appleby, P.G.; Oldfield, F. The calculation of lead-210 dates assuming a constant rate of supply of unsupported ²¹⁰Pb to the sediment. Catena; 1978; 5, pp. 1-8. [DOI: https://dx.doi.org/10.1016/S0341-8162(78)80002-2]

17. Abril, J.M. Radiometric dating of recent sediments: On the performance of ²¹⁰Pb-based CRS chronologies under varying rates of supply. Quat. Geochronol.; 2019; 51, pp. 1-14. [DOI: https://dx.doi.org/10.1016/j.quageo.2018.12.003]

18. Abril, J.M. Multimodal-TERESA, a ²¹⁰Pb-based radiometric dating model for recent sediments under largely varying rates of supply. Quat. Geochronol.; 2020; 55, 101032. [DOI: https://dx.doi.org/10.1016/j.quageo.2019.101032]

19. Caroll, J.; Lerche, I. Sedimentary Processes: Quantification Using Radionuclides; Elsevier: Oxford, UK, 2003.

20. Abril, J.M. A ²¹⁰Pb-based chronological model for recent sediments with random entries of mass and activities: Model development. J. Environ. Radioact.; 2016; 151, pp. 64-74. [DOI: https://dx.doi.org/10.1016/j.jenvrad.2015.09.018] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26421913]

21. Abril, J.M. Why would we use the Sediment Isotope Tomography (SIT) model to establish a ²¹⁰Pb-based chronology in recent-sediment cores?. J. Environ. Radioact.; 2015; 143, pp. 40-46. [DOI: https://dx.doi.org/10.1016/j.jenvrad.2015.02.008] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25725452]

22. Kershaw, P.J.; Woodhead, D.S.; Malcolm, S.J.; Allington, D.J.; Lovett, M.B. A sediment history of Sellafield discharges. J. Environ. Radioact.; 1990; 12, pp. 201-241. [DOI: https://dx.doi.org/10.1016/0265-931X(90)90024-P]

23. Gray, J.; Jones, S.R.; Smith, A.D. Discharges to the environment from the Sellafield site, 1951–1992. J. Radiol. Prot.; 1995; 15, pp. 99-131. [DOI: https://dx.doi.org/10.1088/0952-4746/15/2/001]

24. Camplin, W.C.; Baxter, A.J.; Round, G.D. The radiological impact of discharges of natural radionuclides from a phosphate plant in the United Kingdom. Environ. Int.; 1996; 22, pp. S259-S270. [DOI: https://dx.doi.org/10.1016/S0160-4120(96)00116-X]

25. Othman, I.; Al-Masri, M.S. Impact of phosphate industry on the environment: A case study. Appl. Radiat. Isot.; 2007; 65, pp. 131-141. [DOI: https://dx.doi.org/10.1016/j.apradiso.2006.06.014] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16934479]

26. Appleby, P.G. Chronostratigraphic techniques in recent sediments. Tracking Environmental Change Using Lake Sediments. Basin Analysis, Coring, and Chronological Techniques; Developments in Paleoenvironmental Research Last, W.L.; Smol, J.P. Kluwer: Dordrecht, The Netherlands, 2001; pp. 171-203.

27. Bruel, R.; Sabatier, P. serac: An R package for ShortlivEd RAdionuclide chronology of recent sediment cores. J. Environ. Radioact.; 2020; 225, 106449. [DOI: https://dx.doi.org/10.1016/j.jenvrad.2020.106449]

28. Hillel, D. Soil and Water. Physical Principles and Processes; Academic Press, Inc.: Orlando, FL, USA, 1971.

29. McCartney, M.; Kershaw, P.J.; Allington, D.J. The Behaviour of ²¹⁰Pb and ²²⁶Ra in the Eastern Irish Sea. J. Environ. Radioact.; 1990; 12, pp. 243-265. [DOI: https://dx.doi.org/10.1016/0265-931X(90)90025-Q]

30. Abril, J.M. On the use of ²¹⁰Pb-based records of sedimentation rates and activity concentrations for tracking past environmental changes. J. Environ. Radioact.; 2022; 244, 106823. [DOI: https://dx.doi.org/10.1016/j.jenvrad.2022.106823]

31. Binford, M.W.; Kahl, J.S.; Norton, S.A. Interpretation of ²¹⁰Pb profiles and verification of the CRS dating model in PIRLA project lake sediment cores. J. Paleolimnol.; 1993; 9, pp. 275-296. [DOI: https://dx.doi.org/10.1007/BF00677218]

32. Iurian, A.R.; Millward, G.; Blake, W.; Abril Hernández, J.M. Fine-tuning of ²¹⁰Pb-based methods for dating vegetated saltmarsh sediments. Quat. Geochronol.; 2021; 62, 101153. [DOI: https://dx.doi.org/10.1016/j.quageo.2021.101153]

33. Botwe, B.O.; Abril, J.M.; Schirone, A.; Barsanti, M.; Delbono, I.; Delfanti, R.; Nyarko, E.; Lens, P.N.L. Settling fluxes and sediment accumulation rates by the combined use of sediment traps and sediment cores in Tema Harbour (Ghana). Sci. Total Environ.; 2017; 609, pp. 1114-1125. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.07.139] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28787785]

34. Klubi, E.; Abril, J.M.; Nyarko, E.; Laissaoui, A.; MBenmansour, M. Radioecological assessment and radiometric dating of sediment cores from dynamic sedimentary systems of Pra and Volta estuaries (Ghana) along the Equatorial Atlantic. J. Environ. Radioact.; 2017; 178, pp. 116-126. [DOI: https://dx.doi.org/10.1016/j.jenvrad.2017.08.001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28818644]

35. Bevington, P.A.; Robinson, D.K. Data Reduction and Error Analysis for the Physical Sciences; 3rd ed. McGraw-Hill: New York, NY, USA, 2003.