ESR Dating Ungulate Teeth and Molluscs from the Paleolithic Site Marathousa 1, Megalopolis Basin, Greece

Abstract

At 37°24′ N 22°8′ E, the Megalopolis Basin lies in the central Peloponnese Peninsula, southwestern Greece. In the Megalopolis Basin at ~350 m amsl, the Paleolithic site, Marathousa 1, sits within a palustrine/lacustrine clastic package between Lignite Seams III and II, that both likely correlate with interglacial periods. At Marathousa 1, immediately below Lignite Seam III, lies a clayey-silty sand layer with a horizon rich in molluscs ranging from ~20–40 cm thick. About 0.8–1.3 m below the shell-rich horizon (SRH), lacustrine silty to muddy sands rich in organic matter yielded Paleolithic lithic artefacts associated with Middle Pleistocene fauna, some with cut marks and possible bone knapping, found within palustrine/lacustrine clastic deposits. Since ESR (electron spin resonance) can date teeth and molluscs aged >2 Ma, two bivalve samples, AM66 and AM65, five subsamples from a cervid molar, AT39, and one subsample from another cervid molar, AT68, were independently dated by ESR from Marathousa 1. To calculate the ages, time-averaged cosmic and time- and volumetrically-averaged sedimentary dose rates were calculated using past water depths and sedimentation rates as determined from paleontological and geological criteria. Found in the SRH in Layer UA2, AM66 and AM65 averaged 488 ± 37 ka, which correlates with MIS 13a. Because the bivalves sat stratigraphically above the artefacts and mammalian fossils, their ages constrain the ESR ages for the teeth deposited below. Lying on the unconformity at the base of Layer UA3c with UA4, and its correlative unconformity at the Layer UB4c/UB5 boundary, sat the dated teeth from large mammals. Because the bones in the Palaeoloxodon antiquus skeleton lay in quasi-anatomical association, the likelihood for fossil reworking on the Layer UB3c/UB4 surface is low. Isochron analysis suggests that using a U uptake model with p = 2 provides the most accurate ages for AT39. With p = 2, AT39 dates to 503 ± 13 ka, while AT68 dates to 512 ± 34 ka. Nonetheless, two to three more teeth and molluscs should be dated to confirm these ages, when more samples suitable for ESR dating are found. Both tooth ages correlate well with early MIS 13, an interglacial period with cooler mean global temperatures compared to MIS 11 or 9. Assuming that the archaeological site formed in one event, rather than as a palimpsest, the data suggest that hominins processed elephant and other faunal carcasses along the shores of a shallow lake or marsh in the Megalopolis Basin at 503 ± 12 ka. Between the two horizons dated here, their sedimentation rate averaged 4.8 ± 1.8 to 7.8 ± 2.9 cm/ka.

1. Introduction

In the Pleistocene, Greece likely provided a relatively easy route for hominins migrating from Africa into Europe. Yet, in Greece, the very rare Lower Paleolithic hominid finds have lacked contextual data or derived from a secondary depositional context [1,2]. Recently, however, mammalian fossils, associated with hominin artefacts, have been found in situ in the Megalopolis Basin [3,4,5]. The new ESR (electron spin resonance) dates reported here for teeth found in context at the open-air site, Marathousa 1, in the Megalopolis Basin (Figure 1), provide the first absolute (chronometric) dates for this site, in which artefacts lie near mammalian fossils on the same surface. Mollusc ESR dates are also reported for shells found stratigraphically above the horizon bearing the archaeological materials, providing a capping age for the sedimentary package.

2. Marathousa 1, Megalopolis Basin, Greece

Located at 37°24′ N 22°8′ E at an elevation ~350 m amsl, the Megalopolis Basin lies in the central Peloponnese Peninsula, southwestern Greece (Figure 1). During the Late Miocene-Pliocene, NE-SW tectonic stretching formed a hemi-graben with faults bounding the eastern edge, creating an intermontane depression. In the Pliocene-Pleistocene, a large lake filled the western side of the basin several times, but periodically shallowed to a swamp. The six sedimentary formations in the Megalopolis Basin reach >250 m thick (Figure 1). Organic-rich sediment accumulated mainly along the western margin, while detrital alluvial, fluvial, palustrine and lacustrine facies were deposited to the east. In the late Pleistocene, the current Alfeios River drainage system drained the lake, eroding terraces into the lacustrine sediment [6].

The early-middle Pleistocene Choremi Formation is subdivided into the Marathousa Member, mainly displaying lacustrine or palustrine sediment alternating with lignites, and the Megalopolis Member, with only fluvial sediment. At Marathousa 1, the lacustrine clay, silt, and sand beds host freshwater bivalves, gastropods, and ostracodes (Figure 2). The detrital beds likely formed during dry, cooler conditions (i.e., in glacial or stadial intervals), while the lignite layers probably formed during much wetter, warmer climates (i.e., in interglacial or interstadial periods; [7,8,9]). Since 1969, open-pit mining has found several paleontological localities in the lignite seams [3].

Discovered with associated lithic artefacts during an archaeological survey in 2013 [11,12], stratified bones protruded for ~100 m along a mining profile within the Marathousa Member. Now called Marathousa 1 (MAR1), systematic excavation has continued since 2014 [3,13].

Marathousa 1 sits within the lacustrine sedimentary package between Lignite Seams III and II (Figure 2), correlated to the Middle Pleistocene [7,8,9]. The two excavations, Areas A and B, lie about 60 m apart, but the two show some lateral facies variations between them. Called UA1 in Area A, but UB1 in Area B, Lignite Seam III caps the sequence that includes a mollusc-rich layer, the unit with the artefacts, elephant and other fossils, and other palustrine and lacustrine units. Meanwhile, Lignite Seam II, called UA7 in Area A, but UB10 in Area B, underlies this sedimentary package. Immediately below Lignite Seam III lies a clayey-silty sand layer rich in molluscs (UA2-UB2) ranging from ~20–40 cm thick (Figure 2; [14]). Approximately 0.8–1.3 m below the molluscan horizon sat the zone that yielded the elephant fossils and artefacts. Called UA3 in Area A, but UB4 in Area B, this unit contains lacustrine silty to muddy sands, rich in organic matter averaging ~30 cm thick. In Area A, this unit yielded a nearly complete elephant skeleton in quasi-anatomical position. Most of the mammalian bones and artefacts lay on the unconformity surface between UA3c and UA4 or between UB4c and UB5 (the red line in Figure 2). In Area B, UB5 is a silty sand layer about 25 cm thick that unconformably overlies UB6, but has no equivalent in Area A. UB6 and its correlative UA4 contain a homogeneous, bluish grey, muddy sand with deformation structures. Except for the lithic artefacts, and a few sandstone and limestone clasts, all the layers lacked gravel. Derived from the erosion of previously deposited lacustrine sediment, mud intraclasts averaging <20 mm occurred in several units. High proportions of silty clay suggest that this sediment was deposited in a low-energy system, such as a marsh or a very quiet lake punctuated by high-energy mud flows at the paleolake shores and close to mudflats [14].

From Areas A to B, a lateral facies change occurs, because Area B sat closer to the paleolake shore. In Area B, subaerial floods deposited more and thicker depositional units than in Area A. In both Areas A and B, a major erosional contact divides the sequence between Lignite Seams II and III into two parts (the red line in Figure 2). Overlying Lignite Seam II (UB10 and UA7) in both areas, Units UB6-UB9 and UA4-UA6 contain subaqueously deposited sandy to silty sediment with generally low organic and carbonate concentrations. In Unit UB6 and its correlative UA4, the bluish mud has characteristic load deformation features. In the subaqueously emplaced Units UB5-UB2 and UA3-UA2 generated by subaerial floods, erosional contacts bound each fining-upward package. They all have dark brown, organic-rich, intraclast-rich silty sands at the bottom, overlain by organic-rich massive muds, occasionally interrupted by grey fine sand laminae. Area B had four fining-upward packages, but Area A only saw two deposited. In both areas, the uppermost unit, UA2 and UB2, contains a shell-rich horizon, called the SRH, interbedded within intraclast-rich silty sand. Units UB2 and UA2 have been correlated by their geochemical and sedimentary characteristics [14].

In Area A, Layer UA3c yielded several well-preserved cranial and post-cranial remains from a single Palaeoloxodon antiquus individual. In Area B, more elephantid bones, including another proximal tibia fragment, lay ~60 m from the main elephant bone accumulation. The two left tibiae confirm that the elephant bones derived from at least two individuals. The other vertebrates recovered included teeth, mandibles, and postcranial remains from cervids, bovids, micromammals, turtles, and birds. Since no bones from either Area A or B show any evidence of rounding, bone transport, if any, must have been very limited in space and time, as suggested by the quasi-anatomical orientation of the elephant in Area A. On the bones, minor cracking and laminar flaking suggests very temporally limited bone exposure or geologically rapid burial by a mudflow at the paleolake margin [14,15].

The Marathousa 1 lithic assemblage comprises flakes, flake fragments, cores, core fragments, retouched tools, and debris [12]. In mint condition with no evidence of having been rolled, many flakes and retouched tools bear traces of possible use wear. With no bifaces and large cutting tools, this industry matches well with several other sites, such as Ficoncella, Isernia, and La Polledrara (Italy), that all yielded assemblages that have small tools made on flakes with no handaxes. At Treugol’naya, south Russia, similar flake-tool assemblages have been dated to 583 ± 25 ka (Marine Oxygen Isotope Stage (MIS) 15) until 364 ± 11 ka in MIS 11 [16]. Like Marathousa 1, some of these sites also preserve evidence for elephant or other megafaunal exploitation by early humans [5,17]. In addition to the direct association between the lithics and the fossils at Marathousa 1, especially those with the elephant bones, cut marks preserved on several elephant and other animal bones, and possible bone knapping [15] indicate that hominins exploited these animal carcasses. The small size of the bone fragments suggests that marrow processing occurred at Marathousa 1. The undisturbed sediment and the lithic artefacts’ mint condition, as well as the excellent bone and other organic preservation, all suggest that this area was likely a swamp [3,13,14].

In the Megalopolis Basin, van Vugt et al. [7,8] correlated the sedimentary package between Lignite Seams II and III with MIS 16, ranging from about 670 to 620 ka. They noted, however, that several thinner lignite or organic-rich horizons could correlate with the insolation maxima. They correlated their Layer g (between Lignite Seam II and III) with Insolation Maxima 60, dated to ~650 ka, and placed the Brunhes/Matuyama Boundary in the detrital interval above Lignite Seam I [7,8]. Based on one ESR date at 370 ± 110 ka for the molluscs from the base of their Layer 6 and palynological analyses at a site south of the Marathousa 1 outcrop, Okuda et al. [9] correlated that same package between Lignite Seams II and III with MIS 14 or 12. They correlated their mollusc horizon, which is not likely the stratigraphic equivalent of Layer UA2/UB2, with Insolation Minima 39, but their ESR age would have correlated better with Insolation Minima 34 and thus, would have yielded a more uniform sedimentation rate. They correlated their Layer F (= Layer g) with Insolation Maximum 48–50 at ~507–532 ka. Overall, their section gave a sedimentation rate of 0.21 mm/year [7,8,9]. The high organic content in UA3-UB4 hints that it may correlate with their Layer F. At Marathousa 1, the large faunal assemblage indicates that the horizon with the elephant fossils and artefacts was likely deposited between 0.9 and 0.4 Ma [15], while the microfauna analyses indicate a Middle Pleistocene age [18].

Using the hydroxyapatite (HAP) mineral in vertebrate teeth, ESR can directly date enamel aged ~10 ka to ~2 Ma with 2–8% precision using the enamel’s stable, radiation-sensitive HAP signal. With calcium carbonate minerals in fossils, like calcite and aragonite, ESR can date molluscan shells from ~5 ka to ~2 Ma with 4–10% precision. Thus, ESR dating with two different minerals, namely in molluscs and teeth, acts as two independent dating systems, capable of confirming or constraining each other’s ages. Therefore, two cervid tooth fragments and two mollusc samples were ESR dated in this study (

Table 1. Samples in the Study, Marathousa 1, Greece.

Number		Location			Sample Type
ESR	ESR	Site	Layer/	Depth	Species/	Tooth/
Analysis	Catalogue	Area	Boundary	(cm)	Taxon	Shell Part
AM65	2014.009a1	A	UA2 2	349.94	bivalve	outer shell
AM66	2014.009a3	A	UA2 2	349.94	bivalve	nacreous layer
AT68	2016.001	A	UA3c/UA4 1	350.19	Cervus elaphus	molar
AT39	2014.008	B	UB4c/UB5 1	350.33	cervid	molar

¹ 

The unconformity between UA3c and UA4 in Area A correlates with that between UB4c and UB5 in Area B (see Figure 2).

² 

The Layer UA2 in Area correlates with Layer UB2 in Area B (see Figure 2).

). Although the layer is described as ‘shell-rich’, most of the shells therein occur as small shell fragments embedded in clay, which cannot be easily aggregated into a dating sample.

3. ESR Dating

Since several authors, including Blackwell et al. [19], Rink [20], Schellmann et al. [21], Skinner [22], have discussed the theory and application of ESR dating, only a brief synopsis will be given here. ESR dating methods use radiation-sensitive signals to date hydroxyapatite in vertebrate tooth enamel and calcium carbonate minerals in molluscan shells. Because ESR signals do not zero easily, unlike thermoluminescence (TL), infrared (IRSL) and other optically (OSL) stimulated luminescence signals, ESR signals can be remeasured as often as needed [19,23,24].

Like TL and OSL, ESR dating detects unpaired electrons trapped in mineral crystal lattices. Since their signal height is directly proportional to the number of trapped electrons and to the radiation that the sample has experienced, the sample’s ESR peak height is proportional to 𝖠_Σ, its total accumulated radiation dose. If 𝖠_Σ and the radiation dose rates can be measured, an ESR age can be calculated. While the radiation dose at which the signal saturates sets the maximum ESR dating limit, the minimum spectrometer detection limit sets the minimum ESR dating limit. Low radiation dose rates produce older minimum and maximum dating limits, and vice versa. ESR can date both teeth and molluscs throughout the Pleistocene from 5–10 ka to 2–4 Ma, with 2–8% precision for tooth enamel and 4–10% for molluscs [19,20,23,24,25,26].

The HAP signal in tooth enamel sits at g = 2.0018. Because the HAP peak does not zero with grinding, light exposure, or Earth surface temperatures, its mean ESR signal lifetime, τ, equals ~10¹⁹ years at 25 °C [27]. Q-band ESR shows that recrystallization does not decrease the HAP’s 𝖠_Σ [24,28]. In molluscs, the carbonate peak at g = 2.0007 is usually the most reliable [29], because its τ averages a few million years at 10°C (see Blackwell [30] & references therein; [26]). Ages as old as 2–13 Ma have been reported for molluscs with this peak (e.g., Vichaidid et al. [31], Blackwell et al. [32]). For dating, a radiation-sensitive ESR signal must have been zeroed at the time of interest, and the accumulated dose, 𝖠_Σ, must be less than the saturation dose. If so, ESR ages can be calculated from:

$${\mathsf{A}}_{\mathrm{\Sigma }}={\mathsf{A}}_{\mathrm{int}}+{\mathsf{A}}_{\mathrm{ext}}={\displaystyle \underset{t_0}{\overset{t_1}{\int }}\left(D_{\mathrm{\Sigma }}(t)\right)}\mathrm{d}t={\displaystyle \underset{t_0}{\overset{t_1}{\int }}\left(D_{\mathrm{int}}(t)+D_{\mathrm{sed}}(t)+D_{\cos }(t)\right)}\mathrm{d}t$$

(1)

where

𝖠Σ	=	the total accumulated radiation dose in the sample,
𝖠int	=	the internally derived accumulated dose component,
𝖠ext	=	the externally derived accumulated dose component,
DΣ(t)	=	the total dose rate,
Dint(t)	=	the dose rate from U, its daughters, and other radioisotopes inside the tooth,
Dsed(t)	=	the dose rate from sedimentary U, Th, and K, around the fossil,
Dcos(t)	=	the dose rate from cosmic sources,
t1	=	the sample’s age,
t0	=	today [33].

The added dose method yields 𝖠_Σ, the accumulated dose (see below).

Table A1. Symbols and Abbreviations used in the Paper.

Symbol	Definition
𝖠Σ	the total accumulated radiation dose in the dating sample or subsample
𝖠Σ,i	the total accumulated radiation dose in the ith subsample in an isochron analysis
𝖠int	the internally derived accumulated dose component in the dating sample
𝖠ext	the externally derived accumulated dose component in the dating sample
DΣ(t)	the total dose rate from all sources for the dating sample
Dint(t)	the dose rate from U, its daughters, and other radioisotopes inside the dating subsample
Dint,i(t)	the internal dose rate inside the ith subsample in an isochron analysis
Dext(t)	the dose rate from all sources outside the dating sample
Dsed(t)	the dose rate from sedimentary U, Th, K, and other radioisotopes around the dating sample
Dcos(t)	the dose rate from cosmic sources affecting the dating sample
t1	the sample’s age
t0	today
[Utooth]	the uranium concentration in the whole tooth
[Uen]	the uranium concentration in the enamel for a dated tooth
[Uden]	the uranium concentration in the dentine for a dated tooth
[Umol]	the uranium concentration in a dated molluscan sample
[Used]	the uranium concentration in the sediment around a dating sample
[Thsed]	the thorium concentration in the sediment around a dating sample
[Ksed]	the potassium concentration in the sediment around a dating sample
[Wen]	the water concentration in the enamel in a dated tooth
[Wden]	the water concentration in the dentine in a dated tooth
[Wcem]	the water concentration in the dental cement in a dated tooth
[Wmol]	the water concentration in the in a dated molluscan sample
[Wsed]	the water concentration in the sediment around a dating sample
[Wsed,pal]	the water concentration in the palustrine sediment around a dating sample
[Wsed,all]	the water concentration in the alluvial sediment around a dating sample
${\overline{W}}_{\mathrm{sed}}(t)$	the time-averaged water concentration in the sediment around a dating sample
ρen	the density of the enamel in a dated tooth
ρden	the density of the dentine in a dated tooth
ρcem	the density of the dental cementum in a dated tooth
ρmol	the density of the carbonate shell in a dated molluscan sample
ρsed	the sedimentary density around a dating sample
p	the U uptake rate used in calculating Dint(t)
EU	the early U uptake model used in calculating Dint(t) with p = −1
LU	the linear (continuous) U uptake model used in calculating Dint(t) with p = 0
RU	any recent U uptake model used in calculating Dint(t), often generally used with p = 10
RUx	the recent U uptake model used in calculating Dint(t) with p = x
$D_{\mathrm{ext}}^{\mathrm{BG}}(t)$	the external dose rate derived from the bulk sedimentary geochemcial analysis
$D_{\mathrm{sed},\beta }^{\mathrm{BG}}(t)$	the sedimentary dose rate from β sources derived from bulk sedimentary geochemical analyses
$D_{\mathrm{sed},\mathsf{\gamma }}^{\mathrm{BG}}(t)$	the sedimentary dose rate from γ sources derived from bulk sedimentary geochemical analyses
${\overline{D}}_{\mathrm{int}}(t)$	the time- and volumetrically averaged internal dose rate
${\overline{D}}_{\mathrm{sed}}(t)$	the time-averaged sedimentary dose rate
${\overline{D}}_{\cos }(t)$	the time-averaged cosmic dose rate
${\overline{D}}_{\mathrm{ext}}^{\mathrm{VG}}(t)$	the time- and volumetrically averaged external dose rate from geochemical and cosmic analyses
${\overline{D}}_{\mathrm{ext}}^{\mathrm{iso}}(t)$	the external dose rate derived from the isochron analysis
τ	the mean ESR signal lifetime, a measure of its stability
κα	the α efficiency factor
(234U/238U)0	the initial 234U/238U activity ratio
MIS	Marine (Oyygen) Isotope Stage

lists all the symbols and abbreviations used herein.

3.1. Internal Dose Rates

In both teeth and molluscs, the internal radiation dose rate, D_int(t), mainly derives from U and its daughter isotopes. Modern teeth lack U, but fossil teeth frequently have U concentrations, [U_tooth], exceeding 1 ppm, as in the Marathousa 1 teeth (e.g.,

Table 2. Radioactivity in the Marathousa I Fossils.

Sample		U Concentrations			Internal dose rates, Dint(t)
		Enamel/	Inner	Outer	Enamel/	Inner	Outer
		Shell	Dentine	Dentine	Shell	Dentine	Dentine
		(ppm)	(ppm)	(ppm)	(µGy/y)	(µGy/y)	(µGy/y)
a. AT39, 2014.008, cervid cheek tooth, Area B 1
AT39en1		-	18.41	-
AT39en2		-	16.02	-
AT39en3		-	11.41	-
AT39en1-3		2.08	15.28	21.60	0.230	0.134	0.190
AT39en4		1.50	26.45	15.71	0.171	0.189	0.113
AT39en5		1.53	14.72	17.55	0.176	0.100	0.120
AT39en6		-	14.78	-
AT39en7		-	19.18	-
AT39en6+7		0.79	16.98	19.98	0.089	0.138	0.162
AT39en8		1.90	32.87	17.87	0.218	0.222	0.121
AT39 mean		1.56	19.23	18.54	0.177	0.156	0.141
	±	0.50	7.08	2.28	0.055	0.048	0.033
b. AT68, 2016.001, cervid molar, Area A 1
AT68en1-9		0.07	12.97	15.19	0.008	0.123	0.145
	±	0.02	0.02	0.02	0.003	0.009	0.012
c. 2014.009, bivalves, Area A 2
AM66		0.01	-	-	0.002	-	-
AM65		0.01	-	-	0.002	-	-
mean		0.01	-	-	0.002	-	-
bivalves	±	0.01	-	-	0.002	-	-
Analytical	±	0.01	0.01	0.01	0.002	0.009	0.009
uncertainties 3	-	0.02	0.02	0.02	0.032	0.016	0.033
Detection	~	0.01	0.01	0.01
limit 3	-	0.02	0.02	0.02
Typical water	~	0.02	0.05	0.05
concentrations 3	-	0.02	0.02	0.02

¹ 

D_int(t) calculated assuming a recent uptake model (RU) with p = 2.

² 

D_int(t) calculated assuming a linear uptake model (LU) with p = 0.

³ 

Uncertainties, typical detection limits, and water concentrations depend on the sample’s mass, mineralogy, and tissue type.

), because the dental tissues and bone absorb U post-depositionally. Although teeth absorb substantial U with time, molluscs can also absorb U, but usually only up to 3–5 ppm. To calculate an accurate D_int(t) in any fossil with [U] > 1 ppm, the U uptake rate, p, must either be measured by coupled ESR-²³⁰Th/²³⁴U, where both dating methods measure a unique age. Alternatively, ages can be calculated by assuming an U uptake model: An early uptake model calculation (EU; p = −1) assumes that the sample absorbed all its U soon after death, yielding the youngest possible age for the fossil. For mid Pleistocene teeth, EU model ages always prove to be too young [19]. A linear uptake model age (LU; p = 0) assumes that the sample absorbed U at a constant U uptake rate since its deposition. For many teeth dating to 50–400 ka, LU ages prove to be the most reliable model compared to results from other dating methods, unless the teeth have experienced multiple U uptake events [34]. A recent uptake model age (RU; p > 0, often p = 10) assumes that a fossil absorbed most of its U late in its burial history. In teeth having had multiple U uptake events, RU models with p = 1–5 usually provide more accurate dates, as is often true for teeth older than 400–500 ka [19].

In ESR isochron analyses, a medium-large herbivore tooth, like AT39, which yields ≥4–6 subsamples, can act as its own dosimeter for both its age and its D_int(t), if it has absorbed a range of U concentrations across the tooth. Plotting the subsample’s individual accumulated doses, 𝖠_Σ,i, vs. their individual internal dose rates, D_int,i(t), may yield an isochron. If the plot yields a straight line, the isochron, its slope equals the sample’s age, t₁, and its y-intercept gives the accumulated dose component due to external radiation, 𝖠_ext, and the time-averaged external dose rate, ${\overline{D}}_{\mathrm{ext}}^{\mathrm{iso}}(t)$ [35]. Since the isochron analysis gives a family of lines with different slopes that all converge at 𝖠_ext, the isochron age and ${\overline{D}}_{\mathrm{ext}}^{\mathrm{iso}}(t)$ depend on the U uptake rate, p is accurate, and one selects the p value that gives the ${\overline{D}}_{\mathrm{ext}}^{\mathrm{iso}}(t)$ that best equals ${\overline{D}}_{\mathrm{ext}}^{\mathrm{VG}}(t).$ Thus, by setting ${\overline{D}}_{\mathrm{ext}}^{\mathrm{iso}}(t)$ = ${\overline{D}}_{\mathrm{ext}}^{\mathrm{VG}}(t).$ the isochron method can estimate p [19]. If the plot does not yield a straight line, the isochron analysis can identify a tooth or subsamples within a tooth, whose p has been altered thanks to U remobilization or leaching, and can give clues on the type of U secondary remobilization that the tooth has suffered [34]. In this iterative process, one starts with the parameters from standard age calculations and creates a family of possible isochrons depending on the value of p. Having chosen the value of p whose external dose best matches the one from geochemical analysis, one then recalculates the ages and replots the isochron. Often, this may take three to five iterations.

3.2. External Dose Rates

The external dose rate, D_ext(t), comprises the radiation absorbed from outside the shell or tooth, namely from the sediment, D_sed(t), and from the cosmos, D_cos(t). Although TL dosimetry and γ spectrometry can measure D_ext(t), such measurements rarely occur close enough to the ESR dating samples in the inhomogeneous sediment that often occurs at archaeological sites. In archaeological sites, the actual dose rates can vary by 5–10% in homogeneous sediment, but from 20–500% from spot to spot in the same inhomogeneous layer (e.g., Dibble et al. [36]). If the TL or γ spectrometry testing occurs before excavation begins, datable fossils may lie more than a metre or even farther away from any dosimeters. If the TL or γ dosimetry occurs after the discovery of datable fossils, the sediment around the samples selected for dating will usually occur at some distance from the dosimeters. Ideally, dosimetry should occur before, and as, the sampling for datable material progresses. When the teeth or molluscs are found, however, roughly half of the sediment, and thus, half the measurable D_sed(t), has been lost. Meanwhile, the measurable D_cos(t) has been significantly increased due to the loss of the sedimentary shielding the sample, especially at open air sites. Thus, the TL or γ spectrometry will be inaccurate should one try to measure D_ext(t) by these methods after the samples have been identified. Therefore, geochemical sampling for neutron activation analyses (NAA) must be used to produce more accurate D_sed(t) values by analyzing many samples in inhomogeneous or thinly layered sediment [23].

The sedimentary dose rate, D_sed(t) (e.g.,

Table 3. Sedimentary Radioactivity at Marathousa 1, Megalopolis Basin, Greece.

Sediment	Location		Concentrations			Sedimentary Dose Rates 1
Sample	Area		[U]	[Th]	[K]	$D_{\mathrm{sed},\beta }^{\mathrm{BG}}{(t)}^{2,4}$	$D_{\mathrm{sed},\mathsf{\gamma }}^{\mathrm{BG}}{(t)}^{2,5}$	$D_{\mathrm{sed},\beta }^{\mathrm{BG}}{(t)}^{3,4}$	$D_{\mathrm{sed},\mathsf{\gamma }}^{\mathrm{BG}}{(t)}^{3,5}$
Type	Locus		(ppm)	(ppm)	(wt%)	(mGy/y)	(mGy/y)	(mGy/y)	(mGy/y)
2014.014	B		4.45	7.44	1.10	0.303	1.031	0.206	0.719
bulk sediment		±	0.02	0.16	0.03	0.029	0.065	0.022	0.060
AT39sed1	B		3.92	7.22	1.71	0.368	1.098	0.251	0.766
bulk sediment		±	0.02	0.53	0.06	0.035	0.073	0.028	0.066
AT39sed2 6	B		17.44 6	4.84	0.11	0.514	2.012	0.350	1.405
attached sediment		±	0.02	0.72	0.02	0.046	0.129	0.037	0.119
Mean near AT39	B		8.60	6.51	0.98	0.396	1.382	0.269	0.965
(n = 3)		±	7.66	1.45	0.81	0.238	0.802	0.163	0.562
Mean near AT39 7	B		4.19	7.34	1.41	0.336	1.065	0.229	0.744
(n = 2)		±	0.38	0.16	0.43	0.046	0.047	0.047	0.093
2014.012	A		2.68	4.63	0.46	0.372	0.584	0.254	0.408
bulk sediment	25 cm > shells	±	0.02	0.11	0.01	0.035	0.036	0.029	0.034
2014.010	A		2.95	5.21	0.75	0.516	0.701	0.334	0.489
bulk sediment	shell layer	±	0.02	0.14	0.02	0.048	0.045	0.037	0.041
2014.013	A		3.27	4.92	0.47	0.417	0.659	0.284	0.460
bulk sediment	shell layer	±	0.02	0.12	0.01	0.040	0.042	0.032	0.038
2014.011	A		4.44	5.91	1.02	0.682	0.942	0.465	0.657
bulk sediment	25 cm < shells	±	0.02	0.13	0.03	0.066	0.060	0.053	0.054
Mean near AM66			3.34	5.17	0.68	0.493	0.723	0.335	0.505
(n = 4)		±	0.78	0.55	0.27	0.112	0.111	0.079	0.082
Typical detection limits 8		~	0.01	0.10	0.002	0.002	0.010	0.002	0.008
		–	0.02	0.20	0.004	0.004	0.020	0.003	0.016
NIST 1633B standard			8.79	25.70	1.950	–	–	–	–
		±	0.02	0.70	0.005	–	–	–	–

¹ Abbreviations: $D_{\mathrm{ext},\beta }^{\mathrm{BG}}(t)$ = the external dose rate from β sources, from NAA

                             $D_{\mathrm{ext},\gamma }^{\mathrm{BG}}(t)$ = the external dose rate from γ sources, from NAA

² Calculated with sedimentary water concentration, W_sed =   10.0      ±    5.0    wt%

³ Calculated with sedimentary water concentration, W_sed =    35.0     ±    5.0    wt%

⁴ Calculated with            enamel water concentration, W_en =      2.0     ±    2.0    wt%

                                                                 enamel density, ρ_en =      2.95    ±    0.02    g/cm³

                                                       sedimentary density, ρ_sed =      2.6     ±    0.05    g/cm³

⁵ Calculated using                         cosmic dose rate, D_cos(t) =      0.0     ±    0.0    μGy/y

⁶ Sample likely contaminated with bone or cementum.

⁷  Calculated without AT39sed2.

⁸ Typical NAA detection limits depend on the sample mass and mineralogy.

), derives from sedimentary U, Th, their daughters, and K surrounding the dating sample (e.g., Adamiec & Aitken [37]). Since the enamel or shell dosed by α particles is removed during sample preparation, only β up to ~2–3 mm, and γ up to ~30 cm, away from the sample, are measured. In archaeological sites, most sediment varies substantially from layer to layer and often across a layer, as within the horizon with the elephant carcasses at Marathousa 1. Since each different sedimentary mineral differs in its trace element concentrations, accurate D_sed(t) calculations must volumetrically average the radiation emitted from each mineral component, including bones in the sediment, and each layer within the sediment around each fossil weighted by its volumetric percentage within that γ sphere 30 cm in radius, and within the β sphere ~2–3 mm in radius away from the fossil [38]. In other tests, volumetrically averaged sedimentary dose rates, ${\overline{D}}_{\mathrm{sed}}(t),$ determined from NAA have given comparable dose rates to those from TL or γ dosimetry (e.g., at Grotte de Contrabandiers, Morocco [36]).

Because any covering material attenuates the radiation reaching the sample, the cosmic dose rate, D_cos(t), depends on both the water and sedimentary cover overlaying the sample. If ~100 m of sediment or ~200 m of water covers the sample, the instantaneous cosmic dose rate, D_cos(t), approaches 0 mGy/y. At the modern sediment surface at Marathousa 1 today, D_cos(0) equals 292.9 µGy/y [39] and ~30 m of sediment, but no water, covers the site. For most of its history, however, the site was buried even more deeply, reaching ~100 m by the end of the Middle Pleistocene. Nonetheless, the samples’ D_cos(t) values did not approach 0.0 mGy/y at its initial deposition or during the first few tens of thousands years thereafter. Moreover, the erosional event that occurred after the Middle Pleistocene would have increased D_cos(t) again compared to that at its maximal sedimentary coverage. Therefore, time-averaged cosmic dose rates, ${\overline{D}}_{\cos }(t),$ were calculated for the samples by first finding individual D_cos(t) values [39] and averaging them over time [33].

Reworking a fossil alters both its ${\overline{D}}_{\mathrm{sed}}(t)$ and ${\overline{D}}_{\cos }(t).$ At Marathousa 1, waves and currents in a large paleolake could have reworked shells, as could have floods or storms. To test for reworking in a sedimentary unit, ideally, four to five randomly chosen teeth should be dated. If their 𝖠_Σ’s, ages, or [U_en], or [U_den] differ significantly, reworking is likely [40], but that test has not been possible yet at Marathousa 1, because not enough sample for more than two shell samples and two teeth has yet been found. When more than one layer can be dated, the agreement between the stratigraphic data and the ages can help to demonstrate that reworking has likely not affected the dated samples.

At Marathousa 1, other geological data shows that any reworking was minimal. Several originally articulated bones occurred in almost anatomical association. Although the elephant skeleton was not found in full anatomical position, it was not widely dissociated [14,41]. This suggests that the teeth likely were not reworked much, if at all.

3.3. Sample Preparation

Samples were prepared using standard protocols for a Class 10,000 clean lab. To reduce cross- sample contamination, all glass- and plasticware were soaked in 6 M HCl(aq) and rinsed ≥15 times with doubly distilled, deionized water to remove all Cl⁻ ions [42].

After measuring their mean dentine and enamel thicknesses with a CD-4C digital calliper, both teeth were split into four to eight subsamples with a diamond-tipped Dremel drill, and each subsample was cleaned of attached sediment or dentine, which was removed and saved for NAA. Enamel thicknesses were measured in 30–50 places with a Mitutoyo IP-C112E micrometer, before and after shaving 20 µm off both sides to remove any externally α-doped enamel. After the enamel was powdered to 200–400 mesh (38–76 μm) in an agate mortar and pestle, it was split into 15–16 identical aliquots, each weighing 20.0 ± 0.1 mg [42].

After drilling off any attached sediment and saving it for NAA, the mollusc thicknesses were measured using the Mitutoyo IP-C112E micrometer. To clean them, molluscs were sonicated with 1.0 vol% acetic acid (CH₃COOH) for 1.0 min, rinsed ≥15 times with doubly distilled deionized water, washed once with methanol to chase the water and once with acetone to chase the methanol, and air-dried overnight. After remeasuring the shell thicknesses, the shells were sonicated as above again, dried, and crushed to 100–200 mesh (76–155 μm) in a vice to minimize free radical production. The powders were sonicated as above using 0.1 vol% acetic acid for 15 s. The sample was weighed into 14–16 aliquots, each weighing 70.0 ± 0.1 mg [25].

For the Marathousa 1 fossils, 𝖠_Σ values were calculated by artificially irradiating 15–16 identical aliquots with successively higher, precisely known, added doses, and plotting their ESR signal heights against the known added doses to yield the growth curve whose x-intercept equals 𝖠_Σ. For maximum precision, 𝖠_Σ values were calculated with 15–16 different added doses, the largest of which was ≥10 • 𝖠_Σ (Figure 3; [43,44,45]). All aliquots were irradiated with a ⁶⁰Co γ source with added doses from 0–2560 Grays for the molluscs, and 0–8000 Grays for the teeth at 0.16–0.185 Gy/s. To remove all short-lived interference signals created during the irradiation, all enamel aliquots were annealed at 90 °C for 3.0 days [28], while shell aliquots were annealed at 90 °C for 1.0 h [25].

In the X band, all aliquots were scanned in a JEOL RE1X ESR spectrometer under a 100 kHz field modulation with an amplitude of 0.1 mT at 2.0 mW for the HAP peak with g = 2.0018, and at 5 mW for the aragonite peak at g = 2.0007. Spectra were scanned over 10.0 mT centered at 336.0 mT, using an 8.0 min. sweep time and a 0.1 s time constant. Gains were set to maximize peak heights. Derivative spectra were stored using EWWIN v. 4.5. Peak heights were measured without deconvolution (contra Duval & Grun, [45]), because the two HAP spectral components show identical growth responses with added radiation [28,44]. In the shells, growth curves for four to eight peak heights were built to check for interference in a plateau test.

In this study, to measure the dose rates, D_int(t) and D_sed(t), all associated sediment, and one mollusc, enamel, dentine, and bone per sample, were powdered to ≤500 mesh and analyzed by NAA (e.g., Blackwell et al. [19]; Blackwell [42]) for U, Th, and K. Other enamels, dentines, and molluscs were only analyzed for U. After 60.0 s irradiation and a 10.0 s delay, U was counted for 60.0 s in a neutron counter. Th and K were counted for 20.0 min. in a γ counter. Th was irradiated for 1.0 h and counted after a 7.0 day delay, but K had a 24–30 h delay and 60.0 s irradiation. To ensure accuracy, all results were calibrated against NIST Standard 1633B [42].

3.4. ESR Data Analyses

By plotting the peak heights against the added dose, 𝖠_Σ’s and their errors were calculated using VFit assuming a saturating exponential fit with 1/I² weighting (e.g., Figure 3). Using Rosy v. 1.4.2 or Data-HPS, the ages, dose rates, and their errors were calculated, correcting for attenuation due to β and γ, but not α, back-scattering, water concentration, mineralogy, density, and thickness for the shell, enamel, dentine, and sediment [46]. To calculate the ages, the α/γ efficiency factor, κ_α, was set at 0.15 ± 0.02 for enamel, whereas the molluscan κ_α was set to 0.10 ± 0.02 [47]. Since the site’s geology suggests that little or no Rn loss could have occurred, Rn loss was assumed to be 0.0 ± 0.0 vol%. The initial U activity ratio, (²³⁴U/²³⁸U)₀, was set to 1.2 ± 0.2 [42]. ${\overline{D}}_{\mathrm{sed}}(t)$ was calculated by volumetrically averaging the radiation from each sedimentary unit within the 3 cm or 30 cm sphere of influence around each sample. Today, D_cos(0) = 292.9 μGy/y. Ramped box models (e.g., Deely et al. [33]) were built to calculate individual ${\overline{D}}_{\cos }(t)$ values for each tooth. The time-averaged sedimentary water concentration, ${\overline{W}}_{\mathrm{sed}}(t),$ for the shell horizon within UA2/UB2 was set at 62 ± 2 vol%, while the calculations assumed that ${\overline{W}}_{\mathrm{sed}}(t)$ = 37 ± 5 vol% for the palustrine or lacustrine sediment where the elephant bones were found sitting on the UA3c/UA4-UB4c/UB5 unconformity. Mean ages were calculated by inversely weighting the values by their errors and using Isoplot 3.7 (www.bcg.org/isoplot_etc/isoplot.html).

3.5. Error Analyses

ESR age accuracy depends on the accuracy for the ESR spectrometer’s calibration, the radiation source’s calibration, the sample’s diagenetic state, the sample preparation, and the accuracy of the assumptions made in the age calculations [23]. In interlaboratory calibration tests, Williams College ESR Lab’s dates had no systematic errors (e.g., Barabas et al. [48]; Wieser et al. [49]). Used for the irradiations, the McMaster Nuclear Reactor’s ⁶⁰Co γ source was calibrated daily against a dosimeter that was calibrated annually with the Canadian AECB national ⁶⁰Co γ standard and was calibrated in 2016 against the German ⁶⁰Co γ national standard. Failing to completely remove the outer 20 μm dosed by the external α dose rates from a sample causes its calculated D_int(t) to be underestimated, and its age to be overestimated. Checking the subsample’s thicknesses in >30 spots before and after each side has been shaved off should minimize any from retaining α-doped enamel.

Age uncertainty depends on more than 30 different parameters in the 𝖠_Σ, ${\overline{D}}_{\mathrm{sed}}(t),$ and ${\overline{D}}_{\cos }(t).$ calculations. On average, 𝖠_Σ, ${\overline{D}}_{\mathrm{sed}}(t),$ and ${\overline{D}}_{\cos }(t)$ each add ~1–3% uncertainty to an ESR age. Both calculation programs and the assumptions used in the age calculations include the most recent values and uncertainty estimates for the many parameters (see above).

4. Results

In this dating study at Marathousa 1, five enamel subsamples from AT39, one from another tooth, AT68, and two bivalve samples, AM66 and AM65, were dated (Table 1). AM66 and AM65 occurred in the SRH within Layer UA2 in Area A, which was traced laterally to Area B, where it is called Layer UB2 (Figure 2). In Area B, the SRH sat ~1.3 m above the unconformity between Layers UB4c and UB5. In Area B, draped by sediment from Layer UB4c (Table 1), AT39, a cervid molar, sat on the UB4c/UB5 boundary, was found (Figure 2). From Area B, the UB4c/UB5 unconformity was traced laterally into Area A, where it separates Layers UA3c and UA4. In Area A, AT68, another cervid tooth, was found sitting on the UA3c/UA4 unconformity (Figure 2).

Giving their uncertainties, using neither a double saturating exponential (DSE) nor a linear + saturating exponential (L+ES) fit yielded any significant difference in the calculated 𝖠_Σ compared to that determined with a single saturating exponential (SSE) fit (contra Duval & Grün [44]), probably because annealing removes any short-lived signals produced during the irradiation [27,28]. Moreover, the single saturating fit also produced the smallest uncertainties for 𝖠_Σ for all the fossils (Figure 3). For AM65 and AM66, using the standard single saturating fit produced the same 𝖠_Σ as calculating them with the plateau method (e.g., Schellmann et al. [21]).

4.1. Dating AM66 and AM65

Four sediment samples within 25 cm of AM66 and AM65 were analyzed by NAA. From above, to below, the shell layer, the U concentrations rose from 2.68 ± 0.02 to 4.44 ± 0.02 ppm, while Th rose from 4.63 ± 0.11 to 5.91 ± 14 ppm, and K from 0.46 ± 0.01 to 1.02 ± 0.03 wt.% (Table 3). This variation required that time- and volumetrically averaged sedimentary dose rates, ${\overline{D}}_{\mathrm{sed}}\left(t\right),$ be calculated. The sedimentary dose rates arising from the shells themselves, however, were almost nihil, averaging <0.001 mGy/y for its β contribution, and <0.003 mGy/y for its γ contribution, assuming 40 vol% sedimentary water concentration. In Layer UA2, detailed sedimentological analyses showed that the shells constituted 18.5 wt.% of the bulk sediment, which had 60–65% porosity. When converted into vol%, these data became 6–8 vol% shells, in 33 vol% clastic sediment with 37–42 vol% porosity.

To calculate ${\overline{D}}_{\cos }(t)$, the shells were assumed to have been deposited in water 0–10 m deep, initially at the edge of a marsh or lake, with another 100 m of fluvial and alluvial sediment accumulated continuously above the shells during the Middle Pleistocene, followed by erosion until the overburden reached the current 30 m. When time- and volumetrically averaged, the shells’ ${\overline{D}}_{\cos }\left(t\right)$ equalled 32 ± 6 μGy/y, and their ${\overline{D}}_{\mathrm{sed}}\left(t\right),$ 444 ± 45 μGy/y (Table 3). In the time- and volumetrically averaged calculations for Layer UA2/UB2, the boundaries in each time slice can vary up to 20 ka without significantly altering the shells’ calculated ${\overline{D}}_{\cos }\left(t\right)$ and ${\overline{D}}_{\mathrm{sed}}\left(t\right)$.

Both AM66 and AM65 had U concentrations, [U_mol], measuring 0.01 ± 0.01 ppm (Table 2c). Under 10–50 power, the shells had no significant porosity. While not all leached shells show porosity, its presence would increase the likelihood that the shells had been leached. Only highly weathered molluscan shells usually contain [U_mol] > 2–3 ppm, which can usually be attributed to secondary U uptake or leaching [19,26,32]. Additionally, one would not expect nacreous shell, like AM66, to have high U concentrations [50,51,52]. In molluscs, signal lifetimes vary significantly depending on the peak and species (e.g., see references in Blackwell, Table 2 [30]). Thus, given the lack of similarity in their response to diagenesis and ESR analyses among various mollusc species generally, the fact that the palustrine bivalves, AM66 and AM65, did not respond in the same manner in diagenesis as Tridacna, the giant bivalve from coral reefs [51], should not be surprising.

Encased in palustrine sediment that acted as an aquitard, the likelihood for U leaching from the shells is low. Because the whole Marathousa Member, and hence, the aquitard, stayed unbreached by erosion until mining began, little, if any, of its interstitial water was likely lost after its deposition. Thus, this drastically lowers the possibility for U leaching. The shells’ low U concentrations made an insignificant difference of <1% among the calculated EU, LU, and RU model ages (

Table 4. ESR Ages for the Marathousa 1 Fossils, Greece.

		Accumulated	External Dose Rates		Standard ESR Age 1,2,3,4
Sample		Dose, 𝖠Σ	${\overline{D}}_{\mathrm{sed}}{\left(t\right)}^{1,2,3}$	${\overline{D}}_{\cos }{\left(t\right)}^{1,2}$	EU	LU	RU2	RU8
		(Grays)	(μGy/y)	(μGy/y)	(ka)	(ka)	(ka)	(ka)
a. AT39, 2014.008, cervid tooth, Area B
AT39en1-3y		649.0	680.0	80.0	241.6	378.6	501.0	670.3
	±	22.5	42.0	8.0	14.9	22.0	29.7	43.5
AT39en4		596.7	681.0	80.0	245.8	377.3	491.0	638.5
	±	20.1	42.0	8.0	14.0	21.0	28.7	41.5
AT39en5		582.9	665.0	85.0	268.6	404.1	514.9	650.8
	±	17.8	41.0	9.0	15.7	22.7	30.3	42.7
AT39en6+7		576.0	666.0	84.0	267.8	401.8	512.1	641.1
	±	16.4	41.0	9.0	13.6	20.8	28.1	41.1
AT39en8		648.6	686.0	80.0	240.5	374.3	497.4	666.2
	±	10.8	42.0	9.0	12.5	18.2	26.1	38.6
AT39 mean		621.7	676.0	81.0	253.3	387.0	503.2	654.0
(n = 2)	±	7.1	19.0	4.0	6.3	9.3	12.7	18.5
b. AT68, 2016.001, cervid molar, Area A
AT39en8		567.5	684.0	87.0	303.3	424.4	512.1	643.4
	±	27.2	42.0	9.0	18.3	26.2	34.2	49.4
c. Bivalve shells, Area A
AM66		233.8	444.0	32.0	487.9	488.8	489.9	490.8
2014.009a3	±	5.9	45.0	6.0	53.2	53.6	53.8	54.0
AM65		233.5	444.0	32.0	486.0	488.1	489.1	490.1
2014.009a1	±	3.5	45.0	6.0	52.2	52.6	52.8	53.0
mean		233.6	444.0	32.0	486.9	488.4	489.5	490.4
(n = 2)	±	3.0	45.0	6.0	37.3	37.5	37.6	37.8

¹ Abbreviations: EU = assuming early U uptake, p = −1.

                            LU = assuming linear (continuous) U uptake, p = 0

                          RU2 = assuming recent U uptake, p = 2

                          RU8 = assuming recent U uptake, p = 8

                     ${\overline{D}}_{\mathrm{sed}}(t)$ = the time- and volumetrically averaged sedimentary dose rate

                     ${\overline{D}}_{\cos }(t)$ = the time-averaged cosmic dose rate

² Calculated using

   pallustine time-averaged sedimentary water concentration, $\overline{W}\mathrm{sed}$ = 37.0 ± 5.0  wt%

       alluvial time-averaged sedimentary water concentration, $\overline{W}\mathrm{sed}$ = 10.0 ± 5.0  wt%

³ Calculated using                                              sediment density, ρ_sed = 2.72 ± 0.02   g/cm³

                                                                                      enamel density, ρ_en = 2.95 ± 0.02   g/cm³

                                                                                   dentine density, ρ_den = 2.85 ± 0.02   g/cm³

                                                                     molluscan shell density, ρ_mol = 2.96 ± 0.02   g/cm³

⁴ Calculated using                              enamel α efficiency factor, κ_α,en = 0.15 ± 0.02

                                                         molluscan α efficiency factor, κ_{α mol} = 0.10 ± 0.02

                                                           initial U activity ratio, (²³⁴U/²³⁸U)₀ = 1.20 ± 0.20

                                                               enamel water concentration, W_en = 2.0 ± 1.0    wt%

                                                            dentine water concentration, W_den = 5.0 ± 1.0    wt%

                                                          molluscs water concentration, W_mol = 5.0 ± 1.0    wt%

                                                              radon loss from the tooth, Rn_tooth = 0.0 ± 0.0    vol%

                                                         radon loss from the molluscs, Rn_mol = 0.0 ± 0.0     vol%

c), especially compared to the uncertainties for their ages.

No ESR spectra for AM66 or AM65 showed any evidence of Mn hyperfine splitting (Figure 4). Although some minor organic peaks did occur, they did not interfere with the dating peak. The ESR spectra showed no calcitic peaks. No significant inflexion points occurred in the growth curve (Figure 3), suggesting that the very tiny minor so-called inflexions likely derived from trap competition or lattice defects (cf. [21]). Ages for the two bivalve samples differed by <2 ka, and averaged 487 ± 37 ka (EU) to 490 ± 38 ka (RU2, RU8), and the LU age provides a median age for AM66 at 488 ± 37 ka (Table 4c). This age’s central value correlates with MIS 13b/13a, only a few thousand ages before MIS 12 began. Moreover, these ages do constrain the youngest possible age for the teeth that both underlie the shell horizon. Thus, using the youngest age for the shells, the teeth must predate MIS 11 at the latest, but more likely must predate late MIS 13. Considering their associated uncertainties, this ESR age agreed with the IRSL ages for the sediment package between Lignite Seams II and III [53].

The dates for AM66 and AM65 reported here that correlate with MIS 13, suggest that Lignite Seam III, which sits just above Layer UA2/UB/2, is younger. Thus, Lignite Seam III may correlate with MIS 11, which is younger than the MIS 15/16 boundary age reported by van Vugt et al. [7,8] for the base of Lignite Seam III. The shell dates reported here, however, suggest that Lignite Seam III correlates better with MIS 11 than MIS 13 (contra Okuda et al. [9]) and Lignite II with MIS 15 or earlier. With this correlation, their pollen zone MP3c with low Quercus would correlate with MIS 12, and their thin lignite Layer F, with MIS 13. This correlation also agrees better with the changes in arboreal pollen, especially Quercus, and other changes in MIS 12–16 in Tenaghi Philippon [54].

4.2. Dating the Teeth

In most archaeological sites, dental cementum is lost from the teeth before found or becomes unidentifiable visually. Cementum tends to average <100 µm in thickness [23,24]. In ungulate teeth, depending on the species, enamel plates may have dentine on both sides, cement on one side and dentine on the other, cement covering dentine on one side and dentine on the other, or have dentine on one side only [55,56]. Within a single tooth, enamel plates often have different tissue orientations. More than half of the plates within most cervid or bovid teeth, however, have dentine on both sides. In life, mastication (i.e., wear) can thin or remove the cementum from tooth surfaces. After deposition in sediment, cementum is usually removed, except for sites with the most pristine preservation. In Marathousa 1, simply extracting the teeth from the sediment could have removed any preserved cementum. Whatever the cause, no cementum was recognized in either tooth dated here (Figure 5).

Around the cervid tooth AT39, three sediment samples were analyzed to measure D_sed(t) (Table 3). Given the textural homogeneity of the sediment in the layers, analyzing more than three sediment samples would not have improved the precision of the time-averaged sedimentary dose rate, ${\overline{D}}_{\mathrm{sed}}\left(t\right).$ Since AT39sed2 was scraped off the tooth, dentine, bone, and/or cementum contaminated its sediment, as is evident from its very elevated U and very low K concentrations compared to all the other sediment samples from Marathousa 1. Therefore, AT39sed2 was not used to assess D_sed(t). Given that aquitards bounded the sedimentary package containing the unit with AT39 and AT68, interstitial water surrounding the teeth likely remained intact around the sedimentary grains for most of their depositional history at near or full saturation. Assuming 35 wt.% sedimentary water, D_sed,β(t) for the clastic components averaged 229 ± 47 μGy/y, and D_sed,γ(t), 744 ± 93 μGy/y (Table 3). When volumetrically averaged assuming that AT39’s horizon in UB4, Area B, contained 5 vol% bone, ${\overline{D}}_{\mathrm{sed}}\left(t\right)$ equalled 676 ± 19 μGy/y, while its ${\overline{D}}_{\cos }\left(t\right)$ equalled 81 ± 4 μGy/y (Table 4a).

Initially, AT39 (Figure 5) yielded eight subsamples, but subsamples AT39en1-3 and AT39en6+7 had to be recombined to obtain enough aliquots for the ESR analyses. Even then, AT39en1-3 still did not have enough aliquots to allow it to be dosed highly enough, which led to it being ramped twice (i.e., the subsample’s aliquots were reirradiated, reannealed, and re-ESR analyzed two times; [19]).

AT39’s [U_en] averaged 1.56 ± 0.50 ppm (Table 2a), and its [U_den], 19.23 ± 7.08 ppm U. After AT39en1-3y’s reirradiation, all five AT39’s 𝖠_Σ’s averaged 621.7 ± 7.1 Grays (Table 4a), but the individual values showed considerable variation across the tooth, hinting that an isochron analysis would possibly yield p, the U uptake rate. Given their 2 σ errors, AT39en1-3y at 649.0 ± 22.5 Grays did not differ significantly from the mean. AT39’s LU ages averaged 387.0 ± 9.3 ka, varying from 374.3 ± 18.2 ka to 404.1 ± 22.7 ka (Table 4a). If accurate, these would have correlated with MIS 11.

In AT39, the effect of Rn loss on the ages for AT39en4 was modelled by recalculating the various model ages and their uncertainties with Rn loss from 0–100 vol% (Figure 6). Significance was determined by using a T-test. Had the tooth experienced Rn loss, changing the Rn loss from 0 to 54 vol% would not have produced any significant change in the LU ages, but at 100% Rn loss, the calculated age would have risen to 482 ± 28 ka. In calculating the ages with p = 2 (see below), Rn loss at 0–80 vol% yielded no significant age changes, but 100 vol% Rn loss would have given an age of 602 ± 38 ka. Thus, in all cases, regardless of the model assumed for the age calculation, Rn loss produced older ages. With the aquitards above and below Layer UB4 to limit groundwater movement, however, the likelihood for any substantial Rn loss remained low until mining had breached the aquitards.

Assuming either EU or LU (i.e., −1 ≤ p ≤ 0) for AT39, however, yielded mean ages that were too young compared to the constraint imposed by the bivalve’s ages. This hints that the tooth likely absorbed U with a recent U uptake model, i.e., with a U uptake rate, p > 0. Fortunately, the variability in AT39’s U concentrations and 𝖠_Σ’s suggested that AT39 might yield an isochron (Figure 7). After several iterations through the process, all lines for AT39’s EU, LU, and RU at p = 2, and RU at p = 10, age data sets all converged at 𝖠_ext in the range 395–412 Grays, as did the lines calculated with RU at p = 1, 4, 8, 10, 12, 16, and 20 (not shown on Figure 7 for clarity).

For each uptake model, all the calculated ages fell very close to the isochron lines. Thus, all the lines had significant correlation coefficients at R² = 95%. Assuming LU (i.e., p = 0; the green data on Figure 7), the isochron gave an age of 255 ± 20 ka, while with RU using p = 6 (the pink data, Figure 7) yielded 875 ± 67 ka. Plotting the LU age data (the dotted brown line, Figure 7), as measured from the standard ESR analyses, gives a line with an age of 387 ka and an intercept at 𝖠_ext = 292 Grays. This line fails to converge with any of the model isochron lines. This suggests that the subsamples within the tooth have experienced some secondary U uptake.

By setting the isochron-generated and geochemically-analyzed external dose rates as equal, i.e., ${\overline{D}}_{\mathrm{ext}}^{\mathrm{iso}}(t)$ =${\overline{D}}_{\mathrm{ext}}^{\mathrm{VG}}(t),$ then, one can find p for AT39. Setting ${\overline{D}}_{\mathrm{ext}}^{\mathrm{iso}}(t)$ =${\overline{D}}_{\mathrm{ext}}^{\mathrm{VG}}(t),$ then, gives a line with an intercept of 𝖠_ext = 421 Grays and an age of 387 ka (the dotted orange line, Figure 7). Thus, the isochron analysis shows that a model with p ≅ 2 best approximates AT39’s actual U uptake (the purple line, Figure 7).

Recalculating the standard ESR data for AT39 using an RU model with p = 2 gives an age averaging 503.2 ± 12.7 ka (Table 4b). The Isoplot analyses gave an identical weighted age with 95% (2 σ) confidence. Using Tukey’s biweight calculation gave a weighted average age of 503 ± 15 ka at 95% confidence for AT39, whereas using the robust median method gave a weighted age of ${501}_{-10}^{+14}$ ka with 93.8% confidence. All these calculated mean ages correlate well with MIS 13, a warmer period in the Middle Pleistocene (Figure 8). Nonetheless, a coupled ESR-²³⁰Th/²³⁴U analysis (e.g., McDermott et al. [57]; Grün & McDermott [58]; Jones et al. [59]; Duval et al. [60]) should be run for both AT39’s enamel and dentine to confirm the value for p. Unfortunately, coupled ESR-²³⁰Th/²³⁴U analyses do not always yield reliable results: Successful coupled tooth dates still occur less than half the time on average. Given that teeth act as open systems under every uptake model, except EU, the ²³⁰Th/²³⁴U dating limit does significantly exceed that for ²³⁰Th/²³⁴U closed system dates. Under ²³⁰Th/²³⁴U open systems protocols, the age calculations become more complex. As teeth become older, ages also become less precise, especially when they approach the ²³⁰Th/²³⁴U dating limit [61,62]. Nonetheless, Duval et al. [60] reported coupled ESR-²³⁰Th/²³⁴U ages as old as 1 Ma. With low [U_en] for AT39, ranging from 0.79 to 2.08 ppm (Table 2), getting reliable ²³⁰Th/²³⁴U became more difficult and less precise. Coupled ESR-²³⁰Th/²³⁴U analyses using both enamel and dentine for several subsamples from AT39 failed to give p values. With a very small sample, and [U_en] = 0.07 ppm, ESR-²³⁰Th/²³⁴U analyses for AT68’s enamel will not work. Although coupled ESR-²³⁰Th/²³⁴U analyses could be attempted again, the chance for success remains low with the technology currently available. In sites where both coupled ESR-²³⁰Th/²³⁴U and isochron analyses have been completed, isochron analyses have generated p values that agree well within errors for those found from coupled ESR-²³⁰Th/²³⁴U analyses (e.g., Blackwell et al. [26]).

Since AM66 and AM65 stratigraphically overly AT39, their mean age at 488 ± 37 ka (Figure 8) must set the youngest possible age limit for AT39. Their mean age also does differs significantly from AT39’s EU model age (i.e., p = −1) at the 95% confidence limit (i.e., considering the 2 σ uncertainties; Table 4a). Therefore, AT39 could not possibly have experienced early U uptake, because using an EU model gives AT39’s mean age (at 253 ± 6 ka) that is much too young by > 100 ky! Nor do teeth from other sites with 𝖠_Σ values similar to these usually show U uptake that follows LU models (p = 0), let alone EU models [16,19,23,24,34,63]. At the 95% confidence limit (i.e., considering the 2 σ uncertainties; Table 4), the bivalves’ mean age also significantly exceeds AT39’s LU model mean age (i.e., with p = 0). At the 95% confidence limit, however, the bivalves’ mean age does not differ significantly from AT39’s mean age assuming an RU model with p = 2 at 503 ± 13 ka. This demonstrates that the dated sequence has not been overturned, as do all the other geological, biological, and taphonomic data.

Although the bivalves’ mean age cannot constrain AT39’s potentially oldest mean age, AT39’s standard ESR mean age does converge at an mean age of 654 ± 19 ka with an RU model at p = 8, which is a high p value given the tooth’s 𝖠_Σ values and its U concentrations. The maximum possible age for AT39 ranges from 705 ± 20 ka with p = 10 to 800 ± 25 ka at p = 20, both of which also seem unlikely given the associated fauna found at the site [7,18]. Moreover, finding teeth with p ≥ 10 at other sites is rare [19]. Given the isochron analyses, however, the ages calculated with p = 2 appear more likely to be the more accurate ages. With p = 2, AT39’s mean age correlates well with MIS 13c (Figure 8).

Although AT68 yielded nine subsamples, the tiny pieces had to be recombined into one subsample, AT68en1-9, to get enough aliquots for irradiation and ESR analysis. Its enamel U concentration averaged 0.07 ± 0.02 ppm (Table 2b), and its dentine, 12.97–15.19 ± 0.02 ppm U. Thus, AT68’s enamel U concentrations were considerably smaller than those seen in AT39’s enamel, while AT68’s dentine samples had moderately less U than in some of AT39’s dentine. AT68, however, showed no features to indicate that it had suffered from any U leaching. Unfortunately, its small sample size, low [U_en] and [U_den] precluded it being analyzed by coupled ESR-²³⁰Th/²³⁴U analyses. AT68’s [U_en] and [U_den] also hint that the groundwater salinity at AT68’s deposition site was lower than that seen near AT39. Alternatively, the spot where AT68 sat may have had somewhat lower salinity surface water than where AT39 lay (cf. Blackwell et al. [34]).

With 𝖠_Σ = 567.5 ± 27.2 Grays (Table 4b), AT68en1-9’s LU age at 424.4 ± 26.2 ka does also not differ significantly from some AT39 subsamples. At p = 2, however, AT68en1-9’s age at 512.1 ± 34.2 ka agrees very well with the mean age for AT39 at p = 2. Without an isochron for AT68, however, estimating p for AT68 is an educated guess. Although both teeth sat on the same unconformity and were covered by the same homogeneous sedimentary layer with similar geochemical and sedimentary features [14], the two teeth may not have experienced the same p, especially given their respective U concentrations. Nonetheless, the assumption that both follow roughly the same p does provide a first age approximation for AT68. If AT68 is plotted on the AT39 isochron plot, AT68 does not lie on the same isochron lines for AT39, suggesting that they did experience somewhat different U uptake histories. Given its low U concentration, AT68’s enamel may even have recently lost U. Without coupled ESR-²³⁰Th/²³⁴U ages, which would not have yielded data for the enamel due to its low U concentration, AT68’s p remain uncertain. The constraint imposed by the bivalves’ stratigraphy and ages, however, suggests that AT68 must also be equal or exceed the bivalves’ ages. Assuming p = 2, AT68 dates to 512.1 ± 34.2 ka, which also correlates with MIS 13.

Given that AT68 and AT39 both lay on the same unconformity, albeit separated by ~60 m, the two teeth can be used to calculate a mean average for the whole horizon. Using Isoplot, a weighted average gives an age of 503.0 ± 11.8 ka. Using Tukey’s biweight calculation, the ages average at 505 ± 12 ka, at the 95% (2 σ) confidence limit, whereas using the robust median method gives an age of ${506.6}_{-16.0}^{+8.3}$ ka, with 96.9% confidence. Using feldspar grains, however, Jacobs et al. [52] found IRSL ages correlative with MIS 11–12 for the lacustrine sediment from the package between Lignite Seams II and III. Considering their associated uncertainties, the ESR enamel ages reported here are somewhat older, but agree better with the pollen data.

4.3. Correlations

In Area A, Panagopoulou et al. [3,13] found a partial Palaeoloxodon antiquus skeleton, likely from one individual. About 60 m from the elephant skeleton, more elephantid bones occurred in Area B, including a second left tibia from a different individual. If some bones in Area B should prove to derive from the elephant in Area A, then natural processes or hominids likely caused some minor reworking. Thus far, however, the amount of bone and lithic movement appears to have been minimal. Nonetheless, to test if AT39 and AT68 were reworked, another two to three teeth should be dated for this layer, if suitable teeth can be recovered. The stratigraphy and the anatomical association of the main elephant skeleton, however, indicate that any reworking must have been very limited [15,41]. Finding more teeth that can be sacrificed for ESR dating, given the significance of the faunal and the archaeological finds, will be difficult.

From the stratigraphy, the ages for the cervid teeth AT39 and AT68 must equal or exceed those for AM65 and AM66. That the shell ages do not depend on an U uptake model immediately rules out both EU and LU ages. Sitting on the unconformity at UB4c/UB5, the cervid tooth, AT39, most likely dates to about 503 ± 13 ka and correlates best with MIS 13b (Figure 8). Found on the base of UA3c, which correlates stratigraphically to the base of UB4c, AT68 likely dates to 512 ± 34 ka, and also correlates with early MIS 13 (Figure 8). Since the unconformity at the base of UA3c is the same unconformity at the base of UB4c [14], it dates to 503.0 ± 11.8 ka at the 95% (2 σ) confidence limit, and again correlates with mid MIS 13. Meanwhile, bivalve shells sampled from Layer UA2, AM66 dated to 488 ± 37 ka, also correlating best with MIS 13a. With ~0.8–1.3 m between the base of Layer UA2/UB2 and the unconformity between Layer UB4c/UB5 and UA3c/UA4, the sedimentation rate averages ~4.8 ± 1.8 to 7.8 ± 2.9 cm/ka.

Although MIS 13 (Figure 8) was an interglacial period, it was one of the weakest during the last 800 ka, and similarly, MIS 14 was one of the weakest glacial periods [54,64]. Unlike later Pleistocene interglacials, all the proxy records indicate that mean global temperatures and sea surface temperatures in MIS 13 were not nearly as hot as those in MIS 11c or MIS 5e, or even as warm as those in MIS 1 or 7e. Moreover, in the terrestrial records, the character of the pollen and other climatic proxies are almost indistinguishable within MIS 13–14. For example, at Treugol’naya in the Caucasus Mt., MIS 13 was cooler than during MIS 11 [16], producing mild, relatively dry conditions. At Tenaghi Philippon, MIS 13–14 shows less Quercus overall, coupled with many short-term cyclic expansions and reductions, especially in the latter 20 ka, suggesting greater short-term climatic instability [54], which would correlate with the deposition of AM66 in the shell horizon. At the site southwest of Megalopolis [9], assuming that their Layer F does correlate with MIS 14, the lignite layer is thin with little Quercus, Olea, and Ulmus, hinting at cooler temperatures than during MIS 15 or 11 in the Megalopolis Basin. Except for the warmest times in MIS 13 from ~483–496 ka, the relative temperatures averaged from 3–5 °C cooler than today. Such temperatures are more comparable to those seen in the cooler later parts of MIS 5, 9, or 11 [64]. The ages for AT39 and AT68 suggest that this archaeological site formed just as the temperatures were rising globally from a short-term minimum in MIS 13d. These warm temperatures would have allowed elephants to exploit this environment, where the marsh or lake would have supplied a year-round or almost year-round water supply for the elephants and the other mammalian fauna.

5. Conclusions

This study analysed five subsamples from AT39, including a cervid tooth from Layer UB4c, one subsample from AT68, another cervid tooth from Layer UA3c, and two bivalve subsamples, AM66 and AM65, from shell-rich horizon, the SRH in UA2 from Marathousa 1, Megalopolis Basin, Greece. Independent of the U uptake model assumed for the calculations, a mean ESR age of 488 ± 37 ka for the bivalve shells in the SRH correlates best with MIS 13a. Since the SRH stratigraphically overlies the archaeological site, the mean age for the shells from ARH constrains the minimum age for the archaeological remains in Layer UB4c/UA3c. Hence, the archaeological horizon at Marathousa 1 likely were deposited on the unconformity at the top of Layer UB5/UA4 before MIS 13a. Isochron analysis suggests that AT39 had experienced some secondary U uptake across the whole tooth. That analysis also suggests that a U uptake rate, p = 2, provides the most accurate ESR ages for AT39. Although coupled ESR-²³⁰Th/²³⁴U ages did not yield ages for AT39’s enamel and dentine, while AT68’s [U_en], [U_den], and sample size precluded successful coupled ESR-²³⁰Th/²³⁴U ages to confirm their U uptake rates, p. To check for reworking, at least two to three more teeth should be dated from the unconformity at the base of UA3c/UB4c, but only two teeth from the deposit are available for dating at this time. Assuming p = 2, AT39’s ages averaged 503 ± 13 ka, which correlates with MIS 13b, while AT68 dated to 512 ± 34 ka and also correlates with early MIS 13. Together, the two teeth average to 503.0 ± 11.8 ka. With ~0.8–1.3 m between UA2 and the base of UA3c-UB4c, their sedimentation rate averages ~4.8 ± 1.8 to 7.8 ± 2.9 cm/ka. Assuming that the archaeological site at Marathousa 1 represents one event rather than a palimpsest, at ~503 ± 12 ka, hominins visited the shallow lakeside or marsh in the Megalopolis Basin and hunted or scavenged the skeletons of large mammals there, including elephantids, cervids, and bovids.