^{1}

^{1}

^{1}

^{2}

^{3}

^{1}

^{1}

^{2}

^{3}

This article was submitted to Quaternary Science, Geomorphology and Paleoenvironment, a section of the journal Frontiers in Earth Science

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Shells of bivalve mollusks serve as archives for past climates and ecosystems, and human-environmental interactions as well as life history traits and physiology of the animals. Amongst other proxies, data can be recorded in the shells in the form of element chemical properties. As demonstrated here with measured chemical data (10 elements) from 12 Arctica islandica specimens complemented by numerical simulations, mistakes during sclerochronological data processing can introduce significant bias, adding a further source of error to paleoenvironmental or biological reconstructions. Specifically, signal extraction from noisy LA-ICP-MS (Laser Ablation—Inductively Coupled Plasma—Mass Spectrometry) data generated in line scan mode with circular LA spots requires a weighted rather than an arithmetic moving average. Otherwise, results can be in error by more than 41%. Furthermore, if variations of seasonal shell growth rate remain unconsidered, arithmetic annual averages of intra-annual data will be biased toward the fast-growing season of the year. Actual chemical data differed by between 3.7 and 33.7% from weighted averages. Numerical simulations not only corroborated these findings, but indicated that arithmetic annual means can overestimate or underestimate the actual environmental variable by nearly 40% relative to its seasonal range. The magnitude and direction of the error depends on the timing and rate of both seasonal shell growth and environmental change. With appropriate spatial sampling resolution, weighting can reduce this bias to almost zero. On average, the error reduction attains 80% at a sample depth of 10, 92% when 20 samples were analyzed and nearly 100% when 100 samples were taken from an annual increment. Under some exceptional, though unrealistic circumstances, arithmetic means can be superior to weighted means. To identify the presence of such cases, a numerical simulation is advised based on the shape, amplitude and phase relationships of both curves, i.e., seasonal shell growth and the environmental quantity. To assess the error of the offset induced by arithmetic averaging, Monte Carlo simulations should be employed and seasonal shell growth curves randomly generated based on observed variations.

Bivalve shells contain a wealth of information on past climates (

Parallel to progress in the field of sclerochronology, advances in analytical techniques enabled an increasingly higher sampling resolution. For example, minimum amounts of shell powder required for light stable isotope analysis

In the present paper, we have mathematically assessed both issues. Firstly, we have studied how annual averages computed from arithmetically smoothed LA-ICP-MS line scan data of bivalve shells differ from the equivalent averages filtered with a weighted moving average that considers the circular LA spot geometry. Secondly, we have quantified the bias that is introduced if changes of seasonal shell growth rates remain unconsidered and arithmetic instead of weighted annual means are computed from intra-annual samples. The present study is partly based on published element chemical data (Na, Mg, Mn, Sr, and Ba) of twelve shells of modern

A total of 12 specimens of modern _{2}O_{3} powder, respectively. In case of samples from NE Iceland and the Faroe Islands, one section was immersed in Mutvei’s solution and used for annual increment width measurements in the hinge (data published in

Overview of shells of

Locality | Specimen id | W. d. (m) | Lat/Lon | Age (yr) | Time (yr CE) |
---|---|---|---|---|---|

Faroe Islands | FO10-10-01-V6AR1 | 20 | 62°09′15.00″N, 007°10′25.00″W | 114 | 1895–2006 |

FO10-10-02-V12AR1 | 20 | 62°09′15.00″N, 007°10′25.00″W | 155 | 1855–2005 | |

Isle of Man | IOM0525475R | 40 | 54°07′15.00″N, 004°52′20.40″W | 250 | 1758–2002 |

IOM0505327R | 30 | 54°18′34.80″N, 004°43′14.40″W | 87 | 1919–2003 | |

IOM0505319L | 57 | 54°08′25.80″N, 004°53′58.80″W | 82 | 1923–2002 | |

NE Iceland | ICE120501AL1 | 9 | 66°09′58.92″N, 015°22′58.92″W | 70 | 1946–2010 |

ICE121401AL1 | 10 | 66°11′28.98″N, 015°20′25.44″W | 176 | 1839–2007 | |

ICE120703AL2 | 11 | 66°11′13.68″N, 015°20′48.30″W | 178 | 1837–2006 | |

ICE120505AL2 | 9 | 66°09′58.92″N, 015°22′58.92″W | 80 | 1936–2008 | |

Gulf of Maine | GOM090797R2 | 83 | 44°26′09.83″N, 067°26′18.05″W | 47 | 1963–2008 |

GOM090803R3 | 83 | 44°26′09.83″N, 067°26′18.05″W | 57 | 1953–2007 | |

GOM090829R1 | 83 | 44°26′09.83″N, 067°26′18.05″W | 60 | 1953–2008 |

The remaining polished shell slab (hinge portion) of each specimen was used for ^{2}. During measurement, the laser penetrated approx. 15 µm deep into the sample surface. To remove potential surface contaminants, prior to measurements, sample surfaces were pre-ablated applying a spot size of 100 µm in line scan mode with 3 Hz, 50 μm/s scan speed at an energy density of about 3 J/cm^{2}. This resulted in removal of approx. 0.6 µm of the sample surface. Background intensities were measured for 20 s prior to each line scan.

Accuracy and precision of the analyses were monitored with NIST SRM 610, USGS BCR-2G and USGS MACS-3 (for respective quality control data ^{11}B), sodium (^{23}Na), magnesium (^{25}Mg), aluminum (^{27}Al), potassium (^{39}K), manganese (^{55}Mn), zinc (^{66}Zn), strontium (^{88}Sr), barium (^{137}Ba), and lead (^{208}Pb). Element-specific RSD% (relative standard deviation) values and average detection limits (computed as 3σ_{background} according to ^{43}Ca served as the internal standard for reference materials and shells. Calcium concentrations for the reference materials were taken from GeoReM, for the shells data given in

LA-ICP-MS data obtained in line scan mode show a much higher scatter (= lower signal-to-noise ratio) than results from single spot analysis along the same transect would reveal. To extract signals from such noisy data (

LA-ICP-MS line scan data (here molar Sr/Ca ratios, grey) of the hinge portion of

Schematic illustration of LA-ICP-MS line scan data discretization. Blue circles denote the laser spots on the sample surface. While the laser beam moves forward (here to the right) and ablates shell material, a new data acquisition cycle starts every 1.45 µm (in the setup used in this study). For simplification of the calculations, it is assumed that measurements are instantaneously completed, while in reality, measurements are integrated over some time interval while the laser continues to travel to the right. As the laser spot diameter (h) measures 55 µm and measurements are completed at regular spatial intervals (w) of 1.45 µm, each rectangular stripe of shell (red) measuring 1.45 × 55 µm is measured 38 times (note, illustration not to scale; for simplification, only eight acquisition cycles depicted). Signal extraction of the noisy LA-ICP-MS line scan data should thus be done with a symmetrical moving average integrating 38 data points. However due to the geometry of the laser beam, each rectangular stripe of shell is crossed 38 times (n) by variably sized LA spot portions (pink; note, for illustration purposes

The conversion of the continuous line scan data with overlapping sampling spots into a new chronology with non-overlapping data is known as discretization. Similar techniques as applied herein have previously been used to increase the spatial resolution of LA-ICP-MS data required to obtain a more precise understanding of the rate at which magma ascends (

Conversion of the digitally filtered element chemical data from the distance to the time axis can be done in two ways which were both assessed here. A common practice (also done in

Seasonal growth model (daily resolution) of

In summary, four different variants of annual chronologies were computed and compared to each other: 1) arithmetic annual means from arithmetically smoothed LA data, 2) arithmetic annual means from weight-smoothed LA data, 3) weighted annual means from arithmetically smoothed LA data, and 4) weighted annual means from weight-smoothed LA data. In contrast to arithmetic smoothing, weighted smoothing considers that some shell portions were measured more often than others. In contrast to arithmetic annual means, the calculation of weighted annual means considers the seasonally varying rates of shell growth.

To simulate the differences between arithmetic and weighted annual averages, numerical models were computed in which the sampling resolution (N = 1 to 100 equidistant samples), the seasonal shell growth curve as well as the data recorded in the shell were varied. To provide a tangible example, it was assumed that the bivalve perfectly (i.e., without any kinetic or vital effects) chronicled the sea surface temperature (SST) history during the growing season in its shell (in the form of “proxy-SST”). Aside from assessing differences between arithmetic and weighted proxy-SST averages, this approach also provided the possibility to contrast the shell-derived averages with the actual SST average during the growing season.

Two different model setups were assessed. In the first model, the published seasonal growth model (daily resolution) of

To compute a weighted annual mean from proxy-SST, the relative proportion of time represented by each intra-annual sample served as a weighting factor by which the proxy-SST value from the respective sample was multiplied. The weighted annual average could then be computed as the arithmetic mean of all weighted intra-annual proxy-SST data. Since no actual chemical measurements were used in this model, it was also necessary to simulate the growth rate-related bias contained in the (hypothetical) shell record and determine the proxy-SST value that would have been measured in a given intra-annual sample that formed at daily changing growth rates (evidently, each intra-annual proxy-SST value is likewise biased toward faster growth and differs from the arithmetic mean computed from the simulated environmental SST). This was compensated by weighting the simulated daily environmental SST data by the simulated amount of shell that formed during the respective day. To do so, each simulated daily SST value was multiplied by the corresponding simulated daily increment width. Then, weighted daily temperatures were arithmetically averaged to provide the weighted proxy-SST average of each intra-annual shell sample.

As demonstrated by the representative example depicted in

Annual Sr/Ca chronologies of

Element chemical chronologies of shells from GOM not only were the shortest, but also showed the largest agreement among the four studied regions (average difference: 0.6 vs. 0.9–1.4%, average of the largest difference: 3.7 vs. 5.3–8.8%; ^{2} = 0.35–0.80) between the offset and ontogenetic age only existed for Na, Al, K, Zn, and Pb (

High-resolution proxy data measured in shells can also be combined into annual averages, for example, to identify cycles and trends through time etc. A critical aspect, however, is how these averages are computed, arithmetically or weighted. The latter technique accounts for seasonally varying growth rates and gives more weight to samples from slow-growing portions relative to such from fast-growing portions. A direct comparison of arithmetic and weighted annual means (computed from weight-filtered LA data) revealed more substantial discrepancies than those caused by the different smoothing techniques (

In general, the differences caused by the smoothing techniques and annual averaging methods were larger for shell portions formed during early ontogeny, i.e., in broader increments and when sampling resolution was higher (

The differences between arithmetic and weighted annual means of samples with variable time-averaging were also assessed by numerical simulations (

Numerical simulation of differences between arithmetic and weighted annual means computed from intra-annual proxy-SST.

In other habitats and/or for other species, the seasonal growth and temperature curves can be much different from those presented above. For example, the phase relationship between the curves can vary, the amplitudes can vary, the curves can be compressed or stretched, skewed etc., All of these variations would affect the difference between weighted and arithmetic proxy-SST means as well as their offset from the actual SST mean measured in the environment. Such effects were simulated in an exemplary way and depicted in

Numerical simulation of differences between arithmetic and weighted annual means computed from intra-annual proxy-SST assuming shell growth and SST curves are in phase and have a symmetrical shape.

Numerical simulation of differences between arithmetic and weighted annual means computed from intra-annual proxy-SST assuming maximum shell growth occurs half-way between SST extremes and both curves are symmetrically shaped.

An increasingly shorter growth period (= compressed growth curve) results in larger differences between the actual and reconstructed SST means, unless a weighted mean is computed from a sufficient number of subannual samples (

Numerical simulation of differences between arithmetic and weighted annual means computed from intra-annual proxy-SST assuming main growing season (fast growth) is limited to about 5 months, shell growth and SST curves are in phase and both have symmetrical shape.

Numerical simulation of differences between arithmetic and weighted annual means computed from intra-annual proxy-SST assuming main growing season (fast growth) is limited to about 3 months, shell growth and SST curves are in phase, growth curve is left-skewed and SST curve has a symmetrical shape.

Numerical simulation of differences between arithmetic and weighted annual means computed from intra-annual proxy-SST assuming main growing season (fast growth) is limited to about 5 months, maximum shell growth occurs half-way between SST extremes and both curves are symmetrically shaped.

Numerical simulation of differences between arithmetic and weighted annual means computed from intra-annual proxy-SST assuming main growing season (fast growth) is limited to about 3 months, maximum shell growth occurs half-way between SST extremes, growth curve is left-skewed and SST curve has a symmetrical shape.

Error reduction by computing weighted instead of arithmetic means from intra-annual proxy-SST data of

As demonstrated here, the arithmetic averaging of bivalve shell-derived proxy data representing unequal sample sizes or different amounts of time can challenge paleoenvironmental interpretations. Specifically, signal extraction from noisy LA-ICP-MS data generated in line scan mode with circular LA spots requires a weighted rather than an arithmetic moving average (

According to exemplarily studied shells of

With increasing chronology length, the discrepancies between annual averages computed from arithmetically and circle segment area-based weight-filtered raw element-to-Ca data tended to become larger (

Based on the findings from measured data (

By far more significant than the bias caused by arithmetic data smoothing is the error induced by arithmetic averaging of intra-annual samples with different time-averaging (

To determine how strongly arithmetic means are biased and to obtain the factors by which individual intra-annual samples need to be weighted (^{18}O_{shell} data are aligned in the order of serial sampling until they closely match the shape of a pseudo-δ^{18}O_{shell} curve computed from instrumental water temperature and δ^{18}O_{water} data (^{18}O_{shell} data along the instrumental temperature curve (^{18}O_{water} data—if not available—can also be inferred from salinity using the local freshwater mixing line (e.g.,

Prior to computing weighted annual averages, it should be verified that this mathematical treatment is useful and does not introduce unwanted bias. Specifically, it needs to be explored if the growth rate extremes occur at the inflection points of the SST curve, i.e., half-way between the seasonal extremes of SST (or whatever environmental variable is reconstructed). At such a configuration in conjunction with similarly shaped, symmetrical and smooth (growth and SST) curves, an arithmetic mean can perfectly agree with the actual SST, whereas a weighted average would deviate, in the extreme, by nearly 8% (

As it has become common practice to obtain single samples from annual increments (one sample per annual increment) of bivalve shells (

For various different reasons (e.g., discernability of intra-annual growth patterns and low isotope sampling resolution in narrow increments) it is typically not feasible to determine the timing and rate of seasonal shell growth in all annual increments from which proxy data were obtained, or even in the same specimens. As a practical solution, a seasonal growth model can be constructed from a selection of annual increments (

Seasonal growth curves can differ among specimens and between calendar years, and vary through ontogeny and with increment width. The first and most important thing to test is whether the duration of the growing season changes through ontogeny. If that is the case, a universally applicable seasonal growth model cannot be developed or would require complicated correction factors for different age classes. A quick way to preclude a substantial change in the duration of the growing season is to test whether annual increment width chronologies of specimens with overlapping lifespans can be crossdated. If that is the case and the relative changes in shell growth rate of young and old specimens co-vary, the duration of seasonal shell growth has almost certainly not changed through ontogeny. Because crossdating works well for

In the simulations presented herein, it has been assumed that the shells grow throughout the year, albeit sometimes at an extremely low rate, but never stop growing. If they do stop, the environmental conditions are not recorded, and as a consequence, even the weighted means cannot provide accurate estimates of the actual environmental averages. The assumption of continuous shell growth, however, very likely does not apply at all temporal scales. Many bivalve species (

LA-ICP-MS line scan data should not be arithmetically smoothed if a circular laser beam geometry was used. In such cases, signal extraction is ideally accomplished with a circle segment area-based weighted moving average. This treatment minimizes sample geometry-related errors that can potentially bias paleoenvironmental interpretations. More substantial error is typically introduced if changes in seasonal shell production rates remain unconsidered and arithmetic instead of weighted annual means are computed from intra-annual samples. Shell portions formed during periods of slow growth need to be given greater weight than those formed during fast growth periods. The relative amount of time required to form the respective shell portion from which the sample was obtained serves as the factor by which the associated proxy data needs to be multiplied. As indicated by numerical simulations, arithmetic averages typically overestimate or underestimate the actual environmental variable of interest by up to nearly 40% relative to its range, whereas weighted averages can reduce the error substantially. On average, the error reduction attains 80% at a sample depth of 10, 92% when 20 samples were analyzed and nearly 100% when 100 samples were taken from an annual increment. Under some exceptional, though unrealistic circumstances, arithmetic means can be superior to weighted means. To preclude the presence of such cases, a numerical simulation is advised based on shape, amplitude and phase relationships of the curves of seasonal shell growth and the environmental quantity. To assess the error of the offset induced by arithmetic averaging, Monte Carlo simulations should be employed and seasonal shell growth curves randomly generated.

The problems addressed in the present study are not limited to LA-ICP-MS data, but likewise affect any chemical, physical or structural property recorded in bivalve shells and even other material that formed at variable rate. Therefore, the findings presented here may have wider implications and can also help to improve averages computed from other materials.

The original contributions presented in the study are included in the article/

BS: conceptualization, data curation, formal analysis, investigation, methodology, project administration, supervision, validation, visualization, writing—original draft preparation, writing—review and editing; SM: sample preparation, data curation, methodology, formal analysis, writing—review and editing; RM-K: resources, validation, writing—review and editing; PB: sample preparation, resources, writing—review and editing; UM: resources, writing—review and editing; AW: sample preparation, resources, writing—review and editing; LF: formal analysis, investigation, methodology, writing—review and editing.

This study has been made possible by a German Research Foundation (DFG) grant to BRS (SCHO 793/23).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

We gratefully acknowledge the help of Hilmar A. Holland, Gudrun Thórarinsdóttir, Siggeir Stefánsson, Erlendur Bogason and Sæmundur Einarsson during fieldwork in Iceland as well as Una Matras and Ragnar við Streym (SP/F StreymDiving) for making divers available to collect shells near Faroe Islands. We thank the captain and crew of the RV Prince Madog and of the scalloper Spaven Môr for assistance with the collection of shells from the Isle of Man. We further wish to thank Johanna Boos and Katrin Böhm for their assistance with shell preparation.

The Supplementary Material for this article can be found online at:

^{13}C Suess Effect in the Irish Sea: Results from the First Multi-Centennial Shell-Based marine Master Chronology

^{13}C from Annual Shell Increments of the Ocean Quahog (

^{15}N Values in the Carbonate-Bound Organic Matrix and Periostracum of Bivalve Shells as Environmental Archives

_{3}to Analyze Individual Foraminiferal Shells

^{th}Century in the Southern Barents Sea Derived from Annually Resolved Shell‐Based Records

^{15}N Trends and Multidecadal Variability in Shells of the Bivalve Mollusk,

^{th}Century AD Tephra-Based Radiocarbon Reservoir Ages for North Icelandic Shelf Waters