⁴⁰Ar/³⁹Ar eruption ages of Turkana Basin tuffs: millennial-scale resolution constrains palaeoclimate proxy tuning models and hominin fossil ages

Questions on the origin of our species (Homo sapiens) have captivated humankind for decades. Hominin fossils unearthed following the 1924 discovery of the early hominin Taung skull (Australopithecus africanus) in South Africa (Dart 1925) have afforded remarkable insights into the antiquity (>4.3 Ma) and complexity of human evolution (e.g. Wood and Leakey 2011; Anton et al. 2014). Important hominin fossil sites are now known from South Africa (e.g. Sterkfontein and Malapa), East Africa (e.g. the Olduvai Gorge, the Turkana Basin, Middle Awash and Konso), northern Chad and Eurasia (e.g. Dmanisi).

The Turkana Basin in NW Kenya and SW Ethiopia (Fig. 1) is particularly significant owing to a combination of extensive fossiliferous sedimentary sequences that record much of early human evolution and an abundance of intercalated silicic tuffs (plus key geomagnetic excursions) that provide critical time constraints (e.g. McDougall and Brown 2006, 2008; McDougall et al. 2012). However, most Early Pleistocene–Pliocene tuffs have associated ⁴⁰Ar/³⁹Ar age uncertainties >20 ka (McDougall and Brown 2006, 2008; McDougall et al. 2012) and age estimates for key geomagnetic excursions, such as the base of the Olduvai Subchron (bC2n, 1968–1925 ka; $Δ$ = 43 kyr), are contentious (Channell et al. 2020). These uncertainties inhibit the temporal distinction of closely spaced tuffs (e.g. units within the Chari Member; McDougall and Brown 2006), as well as millennial-scale orbital tuning of sedimentary palaeoclimate cycles (e.g. Joordens et al. 2011, 2013; Lupien et al. 2018, 2020) and the finer scale resolution of the fossil records. The key to resolving these issues is a higher resolution chronological framework for the Turkana Basin tuffs (e.g. Levin 2015).

The order of magnitude improvement in analytical precision afforded by the new generation of multi-collector mass spectrometers should allow the millennial-scale resolution of Turkana Basin tuffs. However, the increased precision also reveals previously unrecognized age heterogeneities, thereby complicating age assignment. In a number of previous studies, feldspars from volcanic tuffs have produced ⁴⁰Ar*/³⁹Ar age distributions with ‘tails’ towards distinctly older (and sometimes younger) ages (e.g. Lo Bello et al. 1987; van den Bogaard et al. 1989; Andersen et al. 2017; Yancey et al. 2018; van Zalinge et al. 2022). The older ages are conventionally attributed to the presence of feldspar xenocrysts, antecrysts and extraneous (i.e. excess or inherited) Ar, although recent studies have invoked low-temperature storage of crystal cargoes and the partial retention of inherited Ar (e.g. Andersen et al. 2017; van Zalinge et al. 2022). Younger feldspar ages are usually ascribed to Ar loss and/or irradiation/analytical artefacts (e.g. Phillips and Matchan 2013; Andersen et al. 2017). The revelation of such dispersed ⁴⁰Ar/³⁹Ar age distributions has necessitated a re-evaluation of how best to extract eruption ages from these datasets (e.g. Schaen et al. 2021 and references cited therein).

This study focused on four palaeoanthropologically important Turkana Basin tuffs, each characterized by distinct ⁴⁰Ar/³⁹Ar feldspar age distributions: the c. 1.9 Ma Kay Behrensmeyer Site (KBS) (= H2 Tuff, Shungura Formation) and the Malbe (= H4 Tuff, Shungura Formation) tuffs and the c. 1.3 Ma Chari (= Tuff L, Shungura Formation) and Gele tuffs in the Koobi Fora Formation (Fig. 1). The ⁴⁰Ar/³⁹Ar results are complemented by elemental data from pumice glass, tuff glass and pumice feldspars, which are used to verify pumice attribution to the host tuffs and provide insights into the complexities of source magma dynamics. We employ a range of statistical approaches to estimate tuff eruption ages before comparing our results with previously published ages. We then utilize the eruption age results to test palaeoclimate proxy interpretations for strata below the KBS Tuff and update age constraints for the associated hominin fossils.

The Turkana Basin formed c. 4.3 myr ago and is dominated by Plio-Pleistocene river, delta and lake sediments, mostly drained from the Ethiopian highlands to the north (Fig. 1). Originally mapped as three separate formations, the sediments of the Koobi Fora (c. 560 m thick), Nachukui (c. 730 m thick) and Shungura (c. 766 m thick) formations are now known to be contemporaneous and are collated within the Omo Group (de Heizelin 1983; Brown and Feibel 1986; Harris et al. 1988). Ash beds in the Omo Group and Konso Formation, Ethiopia are thought to be sourced from volcanic centres within the axial zone of the Main Ethiopian Rift (WoldeGabriel et al. 2005), but the source volcanic centres have yet to be identified. Although some ashfall deposits have been reported, most of the tuffs in the Turkana Basin were formed by fluvial transport and the deposition of ash (e.g. McDougall and Brown 2006, 2008). The generally sharp basal contacts between the tuffs and the underlying sediments, and limited contamination with detrital sedimentary material, indicate a short timespan between eruption and ash deposition in the basin (McDougall and Brown 2006, 2008). Several ash beds, such the KBS Tuff, are regionally extensive, whereas others (e.g. the Gele Tuff) only occur locally (Cerling and Brown 1982; Brown et al. 2006).

A number of Turkana Basin tuffs contain pumice clasts that host feldspar (typically anorthoclase) crystals, which are the primary target for ⁴⁰Ar/³⁹Ar geochronology. Ages are available for c. 40 Omo Group tuffs, ranging from 4.23 ± 0.01 (2$σ$⁠) Ma (unnamed tuff, Apak Member, Nachukui Formation) to 0.750 ± 0.006 (2$σ$⁠) Ma (Silbo tuff, Chari Member, Koobi Fora Formation) (McDougall and Brown 2006, 2008) (Fig. 1).

The first tuff to be identified in the Turkana Basin was the KBS Tuff, which defines the lower boundary of the KBS Member (Brown et al. 2006) (Fig. 1). The type locality of the KBS tuff is archaeological discovery site FxJj1 in Palaeontological Collection Area 105 of the Koobi Fora region. Along the Koobi Fora and Karari ridges, the KBS Tuff occurs in a sequence of lacustrine and fluvial sediments, although airfall ash has been described from Areas 102, 107 and 110 (Brown et al. 2006). Altered KBS tuff outcrops also occur at several localities in the Nachukui Formation, west of Lake Turkana (Brown et al. 2006). The H2 Tuff in the Shungura Formation, Ethiopia is correlated with the KBS tuff on the basis of geochemistry (Cerling et al. 1979) and age data (Drake et al. 1980) (Fig. 1) and occurs as a c. 0.4 m thick unit containing pumice clasts up to 40 cm in diameter (Brown et al. 2006). Although it has been suggested that the KBS (= H2) Tuff is correlated with the Turoha Tuff (TRT) in the Konso region (Katoh et al. 2000), Brown et al. (2006) noted subtle differences in the tuff geochemistry and questioned this correlation.

Initially misidentified as the KBS Tuff (Curtis et al. 1975), the stratigraphically younger Malbe Tuff outcrops in Areas 102, 105, 106 and 112 at Koobi Fora and in the Nachukui Formation (Cerling et al. 1979; Cerling and Brown 1982; McDougall and Brown 2006). The H4 Tuff in the type area of the Shungura Formation is correlated with the Malbe Tuff on the basis of geochemistry (Cerling et al. 1979) and age data (Drake et al. 1980) (Fig. 1). The Malbe Tuff is located c. 15 m above the KBS Tuff in the Koobi Fora (Feibel 1983; McDougall and Brown 2006) and Shungura (Cerling and Brown 1982) formations.

The Chari Tuff defines the base of the Chari Member in the Koobi Fora Formation (Fig. 1) (Brown and Feibel 1985, 1986). The tuff is c. 4.6 m thick in its type section (Area 5), but also occurs across the Ileret region and the Karari Ridge (Areas 131, 133 and 128). The Chari Tuff contains abundant large pumice clasts up to c. 50 cm in diameter. The tuff occurs as a relatively thin layer (c. 10 cm) in the Nachukui Formation (Harris et al. 1988) and has been identified in the Hominin Sites and Paleolakes Drilling Project WTK13 drillcore (Campisano et al. 2017). The Chari Tuff correlates with Tuff L in the Shuguru Formation (Fig. 1) and with the fallout Bright White Tuff (BWT) in the Konso region (Katoh et al. 2000).

The Gele Tuff, named by Feibel et al. (1989), has a limited distribution, outcropping in Area 4 near the Ileret airstrip in the Koobi Fora Formation (Fig. 1). The tuff is c. 0.5 m thick and is considered to be stratigraphically higher than the Chari Tuff on the basis of the available ⁴⁰Ar/³⁹Ar age data (McDougall and Brown 2006).

Previous geochronology results for the KBS (= H2), Malbe (= H4), Chari (= L) and Gele tuffs are summarized in Figure 2 and Table S1. The ages reported in Table S1 include originally published values as well as recalculated ages based on the decay constants of Min et al. (2000) and the fluence monitor ages of Phillips et al. (2022). The ages plotted in Figure 2 and discussed in the following are recalculated values with 2$σ$ uncertainties.

Early attempts to date the KBS Tuff using K–Ar and ⁴⁰Ar/³⁹Ar analyses of bulk feldspar aliquots from entrained pumice clasts produced highly discrepant ages ranging from 2.67 ± 0.54 Ma (Fitch and Miller 1970) to 1.60 ± 0.10 Ma (Curtis et al. 1975) (Fig. 2; Table S1), sparking the so-called ‘KBS Tuff Controversy’ (see Lewin 1987). Subsequent K–Ar (Drake et al. 1980; McDougall et al. 1980; McDougall 1985), ⁴⁰Ar/³⁹Ar (McDougall 1981) and fission track (Gleadow 1980) studies produced consistent ages of c. 1.88 Ma, largely ending the controversy. More recent ⁴⁰Ar/³⁹Ar analyses of single KBS pumice feldspar crystals yielded a precise weighted mean age of 1.873 ± 0.028 Ma (McDougall and Brown 2006), but with mean squared weighted deviate (MSWD) values for individual samples up to 12.1. This raised the possibility of age dispersion beyond analytical uncertainties, which is now clearly evident from multi-collector mass spectrometry analyses. ⁴⁰Ar/³⁹Ar analyses of feldspars from the potentially correlated TRT in the Konso region have a reported weighted mean age of 1.96 ± 0.06 Ma (Katoh et al. 2000), within the uncertainty of these KBS (= H2) results (Table S1).

Age assignment of the Malbe (= H4) Tuff has proved less contentious. Analyses of bulk pumice feldspar aliquots produced consistent K–Ar and ⁴⁰Ar/³⁹Ar ages of c. 1.85 Ma (Drake et al. 1980; McDougall 1985; McDougall and Brown 2006), all within the uncertainty of the reported KBS Tuff ages (Fig. 2). More recent analyses of single feldspar crystals from Malbe pumice clasts produced a reported weighted mean ⁴⁰Ar/³⁹Ar age of 1.850 ± 0.012 Ma (McDougall and Brown 2006), again with elevated MSWD values (up to 15.3) for individual samples.

Reported K–Ar ages for the Chari (= L) Tuff are consistent at 1.39 ± 0.04 Ma and 1.39 ± 0.05 Ma (Drake et al. 1980; McDougall 1985), within the uncertainty of a K–Ar age of 1.34 ± 0.14 Ma obtained from Tuff L samples (Brown and Nash 1976) and a weighted mean ⁴⁰Ar/³⁹Ar age (Chari Tuff samples) of 1.386 ± 0.010 Ma (MSWD values up to 28) (McDougall and Brown 2006). A marginally older ⁴⁰Ar/³⁹Ar age of 1.44 ± 0.04 Ma was obtained for the correlated BWT in the Konso region (Katoh et al. 2000).

The age of the Gele Tuff is constrained by an imprecise K–Ar age of 1.25 ± 0.06 Ma (Brown et al. 1985) and an ⁴⁰Ar/³⁹Ar age of 1.331 ± 0.008 Ma (MSWD values up to 7.4; McDougall and Brown 2006).

A total of 19 feldspar-bearing, in situ pumice clast samples were selected for analysis from KBS/H2, Malbe/H4, Chari/Tuff L and Gele Tuff localities across the Turkana Basin (Fig. 1; Table S2). The samples included six pumice clasts from five KBS/H2 Tuff localities in Areas 105 and 131 of the Koobi Fora Formation and one clast from the H2 Tuff type locality in the Shungura Formation. Five pumice clasts were recovered from two Malbe Tuff outcrops (Koobi Fora Formation) and one H4 Tuff locality (Shungura Formation) (Table S2). Four (c. 20 cm) pumice clasts were selected from three Chari Tuff exposure localities (Koobi Fora Formation), with one pumice clast collected from Tuff L (Shungura Formation). Three pumice samples were recovered from a single exposure of the Gele Tuff in Area 4, Koobi Fora Formation (Fig. 1, Table S2). In addition, a feldspar concentrate from Chari Tuff sample 81-219 (Area 1) was supplied by IM; this sample was analysed previously by McDougall and Brown (2006), yielding an age of 1.378 ± 0.012 Ma (2$σ$⁠). Details of individual pumice samples are summarized in Table S2.

Pumice clasts collected in the field were initially scrubbed with a wire brush to remove any weathered rinds. Each clast was then crushed manually, sieved, washed and dried. Feldspar crystals and glass shards were separated from all KBS/H2, Malbe/H4, Chari/Tuff L and Gele pumice samples for ⁴⁰Ar/³⁹Ar dating and/or electron probe microanalysis (EPMA). Grain mounts of tuff glass (samples ETH86-279, K80-225, K03-0069 and K82-835; Brown et al. 2006; Gathago and Brown 2006) were obtained from the University of Utah collection for reanalysis by EMPA to confirm the affinity of the pumice clasts to the host tuffs (Table S3). The optically clear feldspar cleavage fragments (10–30 grains per sample) selected for ⁴⁰Ar/³⁹Ar analysis were hand-picked under a binocular microscope from either the 0.5–1 mm or 1–2 mm grain size fractions, avoiding grains with obvious inclusions or signs of alteration. To remove undetected adherent glass or secondary minerals, the feldspar grains were treated in an ultrasonic bath with 5% HNO₃ (5 min), 7% HF (10 min), deionized water (10 min) and acetone (3 min).

Feldspar crystals plus tuff and pumice glass fragments were mounted in epoxy resin and analysed at the University of Melbourne using a JEOL JXA-8530F field-emission electron probe microanalyser (Tables S3, S4). The operating conditions included a beam accelerating voltage of 15 kV, a beam current of 5 nA, a beam diameter of 10 $μ$m and counting times of 20 s on peaks and 10 s on two background positions either side of the peak, except for F and Na (10 s/5 s) and Ba (40 s/20 s). Data reduction was achieved using JEOL software with integrated ZAF matrix correction. The standards analysed included wollastonite (Si, Ca), aluminium oxide (Al), magnesium oxide (Mg), hematite (Fe), Mn metal (Mn), titanium oxide (Ti), K-tantalite (K), jadeite (Na), NaF (F), NaCl (Cl), TiZr (Zr) and benitoite (Ba).

Volcanic glass fragments from all the pumice clasts and representative tuff samples were analysed by EPMA to confirm their compositional affinity (Table S3; Fig. 3) because there is evidence of upwards reworking in some instances (e.g. the Morotut Tuff; McDougall and Brown 2006). Tuff glass (and, to a lesser extent, pumice glass) fragments are often variably hydrated post-deposition, leading to relatively low EPMA totals (<96%) (Brown et al. 2006). The EPMA data shown in diagrams are therefore normalized values based on 100% oxide totals. Brown et al. (2006) found that the least mobile elements were Ca, Fe, Ti, Al, Mn, Mg, Cl, Zr and Ba.

Three spots per feldspar grain (core/interior, intermediate and rim/exterior) were analysed for all KBS and Malbe samples, as well as Chari samples TB17-003-P2/P4, whereas two analyses per grain (core/interior and rim/exterior) were conducted on the remaining Chari and Gele samples (Table S4). The presence or absence of compositional zonation was investigated using the Ca and Ba contents (e.g. Iddon et al. 2019). A subset of Chari (= L) grains with relatively elevated core Ca and Ba concentrations were selected for line-scan analyses (Table S4).

In preparation for neutron irradiation, the feldspar grains were packaged into aluminium foil envelopes. To minimize neutron fluence gradients, care was taken to ensure the grains evenly surrounded a smaller (c. 3 mm diameter) aluminium foil envelope containing c. 15 Fish Canyon Tuff sanidine (FCTs) crystals. The packets were then stacked coaxially in silica glass tubes and irradiated in either the CLICIT facility of the Oregon State University TRIGA (OSTR) reactor (UM#75, 10 MWh; UM#77, 20 MWh; UM#79, 20 MWh; UM#85, 4 MWh; UM#89, 5 MWh; UM#94, 6 MWh) or the US Geological Survey's TRIGA reactor (UM#82, 6 MWh).

⁴⁰Ar/³⁹Ar analyses were conducted at the University of Melbourne using a multi-collector Thermo Fisher Scientific ARGUSVI mass spectrometer linked to a stainless-steel gas extraction/purification line and a Photon Machines Fusions 10.6 CO₂ laser system. Details of the ARGUSVI instrument are described in Phillips and Matchan (2013).

We measured ³⁶Ar on a Compact Discrete Dynode detector, with the remaining Ar isotopes measured on Faraday detectors (H1, AX, L1, L2) with low noise 1 × 10¹³ $Ω$ resistors. Following neutron irradiation, the feldspar crystals were loaded into copper sample holders and placed into the stainless-steel sample chamber with a ZnS cover slip. The extraction line was baked at c. 90°C until the extraction line ⁴⁰Ar rate-of-rise levels had decreased to <1 fA min⁻¹.

A number of samples contained relatively small (c. 0.5 mm) feldspar fragments with modest K contents (Fig. 2; Table S3), precluding precise step-heating analyses. Consequently, single-grain, total fusion ⁴⁰Ar/³⁹Ar analyses were carried out on the majority of samples, although incremental heating analyses (two to four steps) were conducted on larger feldspar crystals from KBS Tuff samples 7822-47A, 7822-38 and 77–107y. The main risk of this approach is a reduced capacity to recognize fragments affected by Ar loss (e.g. Andersen et al. 2017; Deino et al. 2019). However, our selection of optically clear feldspar fragments, and the lack of thermal overprints in the region, minimized this risk.

Air aliquots from an automated pipette system were analysed prior to sample analyses to monitor mass discrimination and detector bias. Gas introduced into the ARGUSVI mass spectrometer was equilibrated for 20 s before multi-collector analysis of the five Ar isotopes. Peak signals were collected for a period of 300 s and linearly regressed to time zero (i.e. the time of gas inlet into the mass spectrometer). Extraction line blank behaviour was analogous to that described by Heizler et al. (2021). Mass 40 and 36 blanks stem mostly from the sample chamber and typically decrease after bakeout, with ⁴⁰Ar levels mostly <2 fA (Table S5). ⁴⁰Ar and ³⁶Ar line blanks were measured before every one to three sample analyses, with correction values interpolated from bracketing blank measurements. ³⁷Ar, ³⁸Ar and ³⁹Ar blank levels were minimal, dominated by mass spectrometer contributions and generally time-invariant. These blanks were either averaged across analytical cycles or interpolated where the levels increased or decreased.

Interference correction values for all irradiations, based on analyses of irradiated K-glass and CaF₂ samples in associated (longer) irradiations, are listed in Table S5. Contributions from ³⁶Ar_Cl were calculated using the average of the ³⁶Cl/³⁸Cl production ratios (257.8 ± 2.5) reported for OSTR by Renne et al. (2008) and the (³⁶Ar/³⁸Ar)_Air value (5.3050 ± 0.0084) of Lee et al. (2006). Apparent ages were calculated based on an astronomically tuned FCT sanidine age of 28.176 ± 0.011 Ma (±0.032 Ma, including external uncertainties; Phillips et al. 2022) and the total K-decay constant of 5.463 (± 0.107) × 10⁻¹⁰ a⁻¹ recommended by Min et al. (2000). Unless otherwise stated, age uncertainties are reported at the 2$σ$ level and include uncertainties in the J-values. Full external uncertainties that include uncertainties in the astronomically calibrated age of FCTs and decay constants are reported with the final ages.

The current ⁴⁰Ar/³⁹Ar data were initially filtered to exclude a few analyses with very low Ar yields (<20 fA), low radiogenic ⁴⁰Ar* (<70%), anomalous blanks (>2 fA) or ³⁹Ar signals exceeding detector limits (Table S5). As the tuff feldspars produced variably dispersed ⁴⁰Ar/³⁹Ar data arrays (Table S5), eruption ages were estimated using several statistical approaches trialled in recent studies (see review by Schaen et al. 2021; van Zalinge et al. 2022).

The first approach involves using a ‘gap’ algorithm to eliminate distinctly older ages, followed by data filtering via the normalized median absolute deviation (nMAD) method (e.g. Powell et al. 2020), thereby reducing scatter and enabling the calculation of a weighted mean age (e.g. Deino et al. 2019). This method assumes that the data have a normal (Gaussian) distribution, with older ages attributed to xenocryst/antecryst contamination and younger ages ascribed to Ar loss or irradiation/analytical artefacts.

Schaen et al. (2021) summarized three alternative methods, here designated Models 1–3. Models 1 and 2 assume that the youngest coherent age population most closely approximates the time of eruption. Model 1 (low MSWD weighted mean method) eruption ages are calculated from the youngest group of ages with an MSWD value less than the critical MSWD [= 1 + 2 * SQRT (2/f); f = degrees of freedom]. Model 2 (weighted mean filter) involves calculating the youngest population with a weighted mean age that differs from the next oldest age at the 95% confidence level (e.g. Andersen et al. 2017). Model 3 (normality test and goodness-of-fit filter) assumes that the eruption age is given by the largest subset of data that conforms to a normal distribution, which may not include the youngest ages.

In contrast to these methods, which attempt to extract or define a normal distribution from ⁴⁰Ar/³⁹Ar data arrays, van Zalinge et al. (2022) utilized a Bayesian estimation approach to calculate ⁴⁰Ar/³⁹Ar eruption ages based on the algorithm used for U–Pb zircon age estimations (Keller et al. 2018). They showed that the age arrays obtained from South American ignimbrites followed an exponential distribution, with ⁴⁰Ar/³⁹Ar Bayesian age estimates in general agreement with the U–Pb zircon ages from the same volcanics.

Here, we applied all five approaches to the current data and compared the resultant eruption age estimates. Weighted mean and MSWD values were determined using Isoplot (Ludwig 2012). For the Bayesian age estimates, we used the interactive Jupyter notebook located at https://github.com/brenhinkeller/BayeZirChron.c (Keller et al. 2018), with the data fitted to an exponential distribution (B. Keller, pers. comm. 2020; see worked example in Fig. S1).

The average tuff and pumice glass compositions indicate that the source magmas were peralkaline rhyolites (the KBS/H2 and Chari/L tuffs) and rhyodacites (the Malbe/H4 and Gele tuffs) (Table S3). Compositional data obtained from individual tuffs and pumice glasses support the attribution of the pumice samples to the host KBS (= H2), Malbe (= H4), Chari (= L) or Gele tuffs (Fig. 3; Table S3). The KBS/H2, Chari/L and Gele samples are all characterized by restricted compositional variations (Table S3), with only minor differences between individual pumice samples from the respective tuffs (e.g. slightly elevated TiO₂ in Gele sample TB17-021-P2; Fig. 3). By contrast, the Malbe/H4 Tuff and pumice glass compositions span a much broader compositional range (Fig. 3). Glasses from individual Malbe/H4 pumice clasts exhibit distinct compositional arrays compared with the host tuff, with the latter data spanning the entire compositional range of pumice (Fig. 3). Brown et al. (2006) identified two compositional modes for the Malbe Tuff, whereas the current results indicate a continuum of compositions.

The feldspar compositions for the KBS/H2 Tuff samples define a narrow range (Or_0.34–Or_0.39) and plot close to the anorthoclase–sanidine boundary (Fig. 4). Two feldspar crystals from sample 7822-38 have more sodic compositions, but still plot within the anorthoclase field (Fig. 4a). These two grains also show minor core-to-rim zonation, with the rims exhibiting marginally higher Ca/K ratios. Overall, the restricted range in feldspar compositions mirrors that observed for the KBS/H2 glass compositions.

Feldspar compositions from the Malbe/H4 Tuff pumice samples define an even narrower range (Or_0.34–Or_0.38) and plot adjacent to the anorthoclase–sanidine boundary, mostly within the anorthoclase field (Fig. 4b). However, there are subtle variations between samples, with TB17-022-P1 feldspars exhibiting slightly higher Ca contents and F107-P1 feldspars displaying lower Ca contents. In this case, the limited feldspar compositional range contrasts sharply with the wide range in tuff and pumice glass compositions.

The Chari/Tuff L and Gele feldspars from individual pumice clasts display far more compositional variability than the KBS/H2 and Malbe/H4 samples (Fig. 4c, d). These samples show distinct compositional trends that are near-continuous for most data ≥An_0.07, with a few Gele feldspars plotting at significantly higher anorthite contents (up to An_0.23), extending into the plagioclase field (Fig. 4d). In addition, a number of Gele feldspars plot in the sanidine compositional field.

The majority of feldspar samples from the KBS/H2 and Malbe/H4 samples contain low Ba contents (BaO <0.05 wt%), with no discernible zonation patterns (Table S4). By contrast, the Chari/Tuff L and Gele feldspars are characterized by BaO concentrations ranging from 0.07 to 0.81 wt% and from 0.10 to 0.85 wt%, respectively, with average values of 0.27–0.34 wt% and 0.37–0.40 wt%, respectively (Fig. S2, Table S4). Elevated Ba contents correlate with higher Na₂O/K₂O (Fig. S2) and Ca/K ratios (Table S4). Core–rim analyses of these feldspars reveal some variation within individual fragments, but 73% of Chari/Tuff L grains show no variation in their Ca/K ratios and 61% also show no variation in their Ba contents (Table S4). The remaining Chari/Tuff L grains exhibit discernible core-to-rim zonation, but no consistent zonation patterns – for example, similar proportions of grains show elevated Ca/K and Ba in cores v. rims. By comparison, only 16% of the Gele feldspars are unzoned with respect to their Ca/K values and 44% show no change in their Ba content. Again there is no distinct pattern to the zonation (Table S4). The variations in Ba are highlighted in line-scan analyses completed for three feldspar crystals from Chari sample 81-219 and two feldspar crystals from sample Tuff L-P1, plus a random grain from Tuff L-P1 showing no variation (Fig. S3).

Data from ⁴⁰Ar/³⁹Ar analyses conducted on individual pumice feldspar crystals from the KBS (= H2), Malbe (= H4), Chari (= L) and Gele tuffs are tabulated in Tables S5 and S6 and illustrated in Figures 5 and 6. Estimated eruption ages were calculated using the statistical approaches outlined earlier and are summarized in Table 1. As noted in previous studies of tuff feldspars (e.g. Chen et al. 1996), the data define wedge-shaped patterns in inverse isochron space (Fig. S4).

⁴⁰Ar/³⁹Ar data were obtained from a total of 108 anorthoclase feldspar crystals from six KBS/H2 samples (7822-47A/B, 7822-38, 7722-107, TB17-009-P1, TB17-009-P2 and F102-P) (Tables S5, S6; Fig. 5). The step-heating results from several individual feldspar crystals representing three KBS samples (7822-47A, 7822-38 and 7722-107) show no discernible evidence of Ar loss or gain (Table S5). All the sample datasets exhibit positively skewed age distributions with a ‘tail’ towards older ages and a total range of 1899.1 ± 1.8 ka (2$σ$⁠) to 1872.6 ± 8.4 ka (2$σ$⁠), spanning 36.7 kyr (Fig. 5).

Eruption age estimates for the KBS/H2 samples are tabulated in Table 1. In the case of sample TB17-009-P2, ages determined using the Model 1 and 2 methods are based on only one or two data points. Consequently, the ages were also estimated by excluding the youngest age for this sample. Overall, the nMAD method gives the oldest or equal oldest numerical ages, with the Bayesian estimation approach typically yielding the youngest ages (Table 1; Fig. 5). For five datasets (7822-47B, 7822-38, 7722-107, TB17-009-P2 and F102-P), ages from the various methods are within uncertainty (Table 1). The remaining two datasets show more discordant results, although all age estimates are within 5 ka of one another (Fig. 5).

Weighted mean ages calculated for the seven KBS/H2 samples are reasonably consistent across all methods and range from 1878.9 ± 1.1 ka (2$σ$⁠; MSWD = 2.3; Bayesian method; n = 7) to 1881.2 ± 1.9 ka (2$σ$⁠; MSWD = 5.1; nMAD method; n = 7), a difference of 2.3 kyr (Table 1). The precision of the weighted mean Bayesian age estimate improves to 1879.1 ± 0.6 ka (2$σ$⁠; MSWD = 1.3; P = 0.23) if the youngest age from sample 7822-47A is omitted (Table 1).

As the J-value calculation is a minor contributor to feldspar age uncertainties, we also calculated ages for the combined dataset (Fig. 6), noting that this is not a strictly hierarchical approach to uncertainty estimation (e.g. Renne et al. 1998). The resulting composite age estimates for all methods are indistinguishable from the weighted mean values and are insensitive to the inclusion or exclusion of J-value uncertainties associated with individual analyses. For example, the composite Bayesian age estimate of 1879.1 ± 0.6 ka is identical to the weighted mean value given earlier (Table 1).

The KBS feldspar grains have ⁴⁰Ar/³⁹Ar derived Ca/K values of 0.00762 ± 0.00034 to 0.02068 ± 0.00004 (c. An_0.2–0.7). Notably, there is no obvious correlation between the Ca/K values and apparent ages (Fig. 7).

Single feldspar crystals (n = 76) from five Malbe/H4 Tuff samples yielded ages ranging from 1835.2 ± 2.0 to 1902.5 ± 2.4 ka (2$σ$⁠) (Fig. 5, Tables S5, S6), spanning 71.7 kyr. If the two oldest ages (samples TB17-022-P2 and 81-201) are excluded from consideration, the ‘maximum’ age decreases to 1865.7 ± 1.8 Ma, with a range of 34.3 kyr, analogous to that of the KBS/H2 samples.

Age estimates calculated using the five statistical approaches are summarized in Table 1. Feldspars from two samples (TB17-022-P1 and F107-P1) show concordant results for all methods. For the remaining samples (TB17-022-P2, TB17-022-P3 and 81-201), the nMAD method gives distinctly older age estimates, with Model 3 also yielding older ages for two samples (TB17-022-P3 and 81-201). The remaining methods (Models 1/2, Bayesian) give younger, more consistent ages when considering the uncertainties. Overall, the age estimates vary from 1845.8 ± 1.5 ka (TB17-022-P2; nMAD) to 1835.7 ± 2.1 ka (TB17-022-P3; Bayesian method), a range of 13.7 ka (Table 1).

Weighted mean ages calculated for the five Malbe/H4 samples extend from 1843.1 ± 3.8 ka (MSWD = 19; nMAD method) to 1837.1 ± 0.9 ka (MSWD = 1.4; Bayesian method). The Bayesian approach produces the most concordant age estimates of 1837.4 ± 0.9 ka (MSWD = 2.3; P = 0.06; excluding older outliers) or 1837.7 ± 0.9 ka (MSWD = 1.5; P = 0.19; excluding older outliers and youngest age), with the remaining methods characterized by a low probability of fit parameters (P < 0.01) (Table 1).

As with the KBS/H2 data, we calculated ages for the composite dataset, noting the general concordance of the results with the inter-sample weighted mean determinations (Fig. 6, Table 1). At this juncture, we define a ‘continuum’ group of ages (including individual results within 2.5$σ$ of adjacent data points) because this group shows an improved exponential distribution fit (Fig. 6). For the Malbe/H4 samples, the ‘continuum’ includes all but the three oldest ages. The Bayesian age estimate for the full dataset is 1836.9 ± 1.0 ka (n = 76), regardless of the inclusion/exclusion of J-value uncertainties for individual data points and within uncertainty of the ‘continuum’ group value of 1837.2 ± 0.9 ka (n = 73), indicating the insensitivity of the Bayesian algorithm to older outliers (Table 1).

The Malbe/H4 feldspar grains have ⁴⁰Ar/³⁹Ar derived Ca/K values of 0.00517 ± 0.00034 to 0.07089 ± 0.00034 (c. An_0.2–2.1). In accord with the KBS/H2 results, there is no obvious correlation between the Malbe/H4 Ca/K values and apparent ages (Fig. 7).

Chari Tuff and Tuff L feldspars from five samples show strongly positively skewed apparent age distributions (Figs 5, 6). If all the data are aggregated (n = 114), the ages range from 1347.3 ± 5.6 to 1949.4 ± 4.0 ka (2$σ$⁠), an overall timespan of 611.7 kyr (Table S5). If ages older than the ‘continuum’ group are excluded (n = 17), then the ‘maximum’ age becomes 1436.9 ± 1.6 ka, a span of 85.6 kyr (Table S5).

Quantile–quantile plots (Fig. S4) and other tests of normality show that the data in the Chari/Tuff L age continua (i.e. excluding older outliers) are not normally distributed, such that the calculation of the weighted mean or median absolute deviation values is inappropriate. This is clearly evident from Table 1, where the calculation of nMAD and Model 1–3 ages requires the exclusion of a majority of data points, yielding discordant weighted mean ages with low probabilities of fit (<0.01). By contrast, the Bayesian age estimates for the six datasets (all data points) have a more limited dispersion, ranging from 1351.3 ± 6.3 ka (K82-787b) to 1364.7 ± 4.2 ka (81-219) (Table 1). The continuum datasets show slightly more precise age estimates of 1351.9 ± 6.3 ka (K82-787b) to 1364.1 ± 2.8 ka (81-219), for a weighted mean age of 1359.3 ± 3.8 ka (MSWD = 3.8; P = 0.00). If the oldest age estimate from sample 81-219 is omitted, a more precise weighted mean age of 1357.5 ± 1.8 ka (MSWD = 0.9; P = 0.45) is calculated.

The Bayesian age estimate for the combined dataset is 1356.5 ± 1.7 ka (n = 113; excluding the youngest age from K82-787b) and 1357.1 ± 1.7 ka for the continuum age array (n = 96) (Table 1). As with the KBS/H2 and Malbe/H4 samples, these values are indistinguishable from the inter-sample weighted mean age estimates calculated earlier (Table 1).

The majority of the age continuum grains from the Chari/Tuff L samples have Ca/K values <0.05 (c. An_1.7), with values for the remaining grains ranging up to 0.18 (c. An₅; Fig. 7). Feldspar fragments forming the older ‘tail’ of the age distribution have Ca/K values of 0.014 ± 0.002 to 0.873 ± 0.004 (c. An_0.5–13.2) and there is a discernible positive correlation between the Ca/K values and apparent ages (Fig. 7, Fig. S4).

Feldspars from the Gele pumice clasts show a similar positively skewed age distribution to the Chari/L samples. Overall, the apparent ages (n = 48) range from 1306.9 ± 5.2 ka (sample TB17-021-P1) to 1570 ± 11 ka (sample TB17-021-P2), a range of 279.3 kyr (Table S5; Figs 5, 6). The continuum age array (n = 41) has a maximum age of 1389.2 ± 6.0 ka (sample TB17-021-P2), giving a reduced age difference of 93.5 kyr.

As per the Chari/L data, the Gele age continuum data are not normally distributed (Fig. S5), thus complicating the calculation of the weighted mean or median absolute deviation age estimates. This is again evident from the small number of data points included in the age estimates and the low probabilities of fit (<0.01) for most methods (Table 1).

Bayesian age estimates for the three Gele samples vary from 1309.0 ± 6.0 ka (TB17-021-P1) to 1316.8 ± 9.9 ka (TB17-021-P3). Consideration of the ‘continuum’ data points only gives age estimates of 1310.0 ± 5.6 ka to 1317.6 ± 8.2 ka (Table 1). Weighted mean ages for the three samples are 1314.9 ± 2.3 ka (MSWD = 1.9; P = 0.15; n = 3) for the continuum array and 1315.4 ± 1.9 ka (MSWD = 0.3; P = 0.75) if the youngest age from sample TB-021-P1 is excluded (Table 1).

The composite dataset has a Bayesian age estimate of 1315.4 ± 2.4 ka (n = 47) or 1315.4 ± 2.0 ka (n = 40) for the continuum age group, excluding the youngest data point (Fig. 6). As noted for previous tuffs, these values are unaffected by the inclusion or exclusion of J-value uncertainties and are indistinguishable from the weighted mean age estimates.

The Gele continuum feldspars have Ca/K values of 0.0476 ± 0.0006 to 0.483 ± 0.004 (c. An_1.7–11.0), whereas the grains in the older ‘tail’ have Ca/K values of 0.165 ± 0.001 to 1.208 ± 0.002 (c. An_4.7–15.0) (Fig. 7). As observed for the Chari/Tuff L samples, there is a discernible correlation between the Ca/K values and the apparent ages (Fig. 7).

Feldspar crystals from the KBS/H2, Malbe/H4, Chari/Tuff L and Gele Tuff samples all show ⁴⁰Ar/³⁹Ar age variability beyond that attributable to analytical or reactor-related uncertainties (Figs 5, 6). A similar dispersion has been observed in other silicic volcanic rocks, where feldspar ⁴⁰Ar/³⁹Ar age arrays typically define broadly exponential distributions (e.g. Andersen et al. 2017; van Zalinge et al. 2022). The older ages in these arrays are usually attributed to excess Ar contamination and/or the presence of antecrysts/xenocrysts (Schaen et al. 2021 and references cited therein). By contrast, younger ages are typically considered to be caused by Ar loss, irradiation artefacts (e.g. fluence gradients, self-shielding) or analytical bias (e.g. blank aberrations, hydrocarbon interferences) (Schaen et al. 2021 and references cited therein).

Of the four Turkana Basin tuffs, the KBS/H2 and Malbe/H4 Tuff samples are characterized by broadly similar ⁴⁰Ar/³⁹Ar age distributions with restricted age ranges (37–78 kyr, respectively), narrow feldspar compositional fields (Fig. 4) and no obvious correlation between the feldspar ages and composition (e.g. Ca/K ratios; Fig. 7). By contrast, the Chari and Gele feldspar samples embody age ranges >250 ka (Figs 5, 6), broader feldspar compositional variations (Fig. 4) and distinct correlations between age and Ca/K ratios (Fig. 7). The latter correlations suggest that the older ⁴⁰Ar/³⁹Ar ages are more likely to be due to the retention of inherited (i.e. pre-eruption) radiogenic Ar rather than excess Ar contamination because the latter is unlikely to be related to the compositional trends. This interpretation is supported by recent studies of silicic volcanic systems showing consistent U-series/U–Pb zircon and ⁴⁰Ar/³⁹Ar feldspar ages, indicating the ‘cold storage’ of crystal accumulations for periods up to several million years (Cooper and Kent 2014; van Zalinge et al. 2022). Modelling of radiogenic Ar accumulation and diffusional loss rates in crystal reservoirs indicates that feldspar crystals must be stored below c. 475°C to retain inherited Ar (Andersen et al. 2017; van Zalinge et al. 2022). In addition, magma residence times prior to eruption must be short (years to tens of years) to preserve older ⁴⁰Ar/³⁹Ar ages during exposure to elevated liquidus temperatures (>750°C) (van Zalinge et al. 2022). The less than perfect relationship between age and composition in the Turkana tuff samples is probably a function of local magma temperatures, the time of exposure of crystals to the magma prior to eruption and the diffusivity characteristics of individual feldspar crystals, although we cannot totally discount the possible presence of excess Ar (e.g. from glass inclusions).

This scenario accords with our current understanding of the dynamics of large silicic magmatic systems (e.g. Cooper and Kent 2014; Cashman et al. 2017; van Zalinge et al. 2022), which are viewed as a series of crustal-scale reservoirs dominated by crystal mush and crystal accumulations, stored at temperatures ranging from liquidus values to below c. 450°C (Cooper and Kent 2014; van Zalinge et al. 2022). The limited volumes of melt associated with these reservoirs also accounts for the lack of geophysical evidence for large melt accumulations associated with silicic systems (e.g. Hübert et al. 2018). Intrusions of mafic magmas are thought to produce periodic silicic magma hot zones in the upper crust that occasionally erupt explosively (Cashman et al. 2017; Sparks et al. 2019).

The current combination of ⁴⁰Ar/³⁹Ar and glass/feldspar compositional data provides new insights into the eruption dynamics of the source volcanic systems that produced the Turkana tuffs. The restricted age and chemistry of the KBS/H2 samples (Figs 4, 6) implies the eruption of a homogenous, more evolved rhyolitic magma with very limited inclusions of older antecrysts. The variable glass chemistry of the Malbe/H4 samples (Figs 4, 6) suggests a more heterogeneous evolved rhyodacite magma, possibly caused by the mixing of magmas with slightly different compositions. The tight clustering of feldspar compositions is again consistent with the limited entrainment of older antecrysts (or xenocrysts).

The elevated Ba contents of the Chari/Tuff L and Gele feldspars (Table S4) indicates that their source magmas were less evolved than those that produced the KBS/H2 and Malbe/H4 tuffs because Ba partitions preferentially into early feldspar crystals (e.g. Iddon et al. 2019). The distinct feldspar major element arrays defined by the Chari/L and Gele feldspars (Fig. 4) imply continuous fractional crystallization trends. However, the higher, more variable Ba contents associated with the sodic feldspars (Figs S1, S2) indicate a more complex scenario, with these grains likely to be antecrysts. This accords with similar feldspar chemistries observed in younger peralkaline rhyolitic volcanics from the Main Ethiopian Rift (Iddon et al. 2019). The Chari/Tuff L and Gele eruptions appear to have been generated from relatively homogenous melts that entrained juvenile as well as older, antecrystic feldspar crystals. Broadly similar eruption models have been proposed for the Bishop Tuff (Andersen et al. 2017) and Central Andes volcanics (van Zalinge et al. 2022), suggesting that such eruption dynamics are more typical of silicic magma systems than previously envisaged.

It is generally accepted that tuff eruption ages are best represented by the youngest subset of concordant feldspar ages that are unaffected by inherited Ar, excess Ar, Ar loss or other irradiation/analytical biases leading to anomalous ages (Schaen et al. 2021 and references cited therein). However, there is limited agreement on the best approach(es) for constraining these ages. Differences in eruption age estimates from the five statistical methods trialled in this study are relatively minor (c. 2 ka) for the KBS (= H2) Tuff, which is characterized by the smallest age range and a near-normal age distribution, but are more significant (6–9 ka) for the Malbe/H4, Chari/Tuff L and Gele feldspar samples, which exhibit non-normal distributions (Table 1).

Of the statistical approaches employed, the nMAD method typically produced the oldest age estimates, followed by the normality test and goodness-of-fit filter method (Model 3) (Table 1, Fig. 8). This is not entirely unexpected because the nMAD/gap outlier exclusion approach essentially forces a normal distribution by excluding both older and younger ages (e.g. Deino et al. 2019). Given the distinctly non-gaussian distributions of the current datasets (in particular the Chari/Tuff L and Gele samples), this method typically yields the most discordant intra-sample age results and is thus non-ideal for such dispersed age distributions. The normality test and goodness-of-fit filter method (Model 3) also yields older eruption age estimates in most cases because it is based on the largest array of data points that define a normal distribution and meet goodness-of-fit parameters (e.g. P > 0.05; see Schaen et al. 2021). This method produced the least consistent intra-sample results and may be better suited to larger datasets. The low MSWD weighted mean (Model 1) and weighted mean filter (Model 2) methods invariably give younger age estimates that should equate with eruption events if the data are minimally affected by Ar loss or other irradiation/analytical artefacts that would cause anomalously young ages. However, these methods again show poor inter-sample age consistency for the Turkana tuffs, with younger ages tending to be associated with larger sample sizes (Schaen et al. 2021).

In contrast with these approaches, the Bayesian estimation approach minimizes the need for older outlier identification, removing some subjectivity. The insensitivity of the approach to older outliers is illustrated by the equivalence of ages determined for all data points compared with those calculated from the continuum datasets (Table 1). However, the Bayesian algorithm is associated with the youngest eruption age estimates for most samples and shows some sensitivity to young outliers. Nonetheless, the Bayesian method shows the most consistent inter-sample ages per locality, with concordance further enhanced for the KBS/H2 and Gele Tuff samples by omitting just two anomalously young data points (KBS sample 7822-47A and Gele sample TB17-021-P1) (Table 1). The exclusion of additional young data points has a minimal impact on the ages, suggesting that the two (anomalously young) KBS/Gele grains may be affected by Ar loss and/or irradiation/analytical aberrations.

Despite the advantages of a stochastic age estimation approach, the possibility of age bias towards younger ages due to issues related to Ar loss cannot be ignored. We attempted to mitigate this issue through the careful selection of optically clear feldspar fragments devoid of obvious alteration, combined with evidence from EPMA data and the available (albeit limited) step-heating results, which show no evidence of alteration or Ar loss, respectively. In addition, the concordance of inter-sample Bayesian age estimates suggests that Ar loss is not a major factor. Nevertheless, further work is required to test the Bayesian approach more rigorously for age robustness.

As noted earlier, similar eruption age estimates were obtained for the KBS/H2 feldspar samples from all statistical approaches. However, as the Bayesian age estimates are the most concordant across all samples, our preferred age for the KBS (= H2) Tuff is 1878.9 ± 1.1 ka (MSWD = 2.3; P = 0.03) or 1879.1 ± 0.6 ka (MSWD = 0.61; P = 0.31) if the youngest age (sample 7822-47A) is excluded from consideration. The inclusion of external uncertainties increases these values to ±2.6 ka and ±2.4 ka, respectively (Table 1). Our preferred age for the KBS (= H2) Tuff (1879.1 ± 0.6/2.4 ka) is within the uncertainty of previous bulk K–Ar, ⁴⁰Ar/³⁹Ar and fission track ages (Table 1), apart from the considerably older values reported by Fitch and Miller (1970, 1976), Hurford et al. (1976) and Fitch et al. (1996) (Fig. 2, Table S1). This result is also indistinguishable from the weighted mean ⁴⁰Ar/³⁹Ar age of 1873 ± 28 ka obtained from single-feldspar analyses by McDougall and Brown (2006), despite the use of different statistical approaches. In this case, the similarity is attributable to the restricted age range of the KBS/H2 dataset (Fig. 6).

The ⁴⁰Ar/³⁹Ar age of the TRT in the Konso Formation (1960 ± 60 ka; Katoh et al. 2000) is within the uncertainty of the KBS age reported by McDougall and Brown (2006), but is distinctly older than the current eruption age estimate (1879.1 ± 0.6 ka). Nonetheless, the youngest individual TRT feldspar ages (e.g. 1878 ± 80 ka; 2$σ$⁠) are within the uncertainty, raising the possibility that, if equivalent, the TRT may contain a higher proportion of older feldspar antecrysts. Although Brown et al. (2006) questioned the correlation with the KBS Tuff based on differences in the glass geochemistry, the age and geochemistry of the TRT are not sufficiently dissimilar to rule out a correlation with the KBS/H2 Tuff and further geochronology and geochemistry are warranted to test this relationship.

The Bayesian approach again gives the most concordant age estimates for the Malbe/H4 samples and our preferred ages are 1837.4 ± 0.9 (±2.4 ka; MSWD = 2.3; P = 0.06) or 1837.7 ± 0.9 (±2.4 ka; MSWD = 1.5; P = 0.19), excluding the youngest age (full external uncertainties are shown in parentheses) (Table 1). These values are within the uncertainty of the K–Ar (1850 ± 40 ka) and ⁴⁰Ar/³⁹Ar (1850 ± 12 ka) ages reported by McDougall (1985) and McDougall and Brown (2006), respectively, but slightly older than earlier K–Ar ages of Drake et al. (1980) and McDougall et al. (1980);(Table S1).

Our study reveals that the eruptions of the KBS and Malbe tuffs were separated by c. 40–43 kyr, where the previous chronology was indistinguishable (McDougall and Brown 2006) (Fig. 6, Table S1). The Malbe/H4 Tuff overlies the KBS/H2 Tuff by c. 16 m in the Shungura Formation (Cerling and Brown 1982) and c. 15 m in the Koobi Fora Formation (Feibel 1983; McDougall and Brown 2006), indicating local sedimentation rates of c. 39.7 and c. 37.3 cm ka⁻¹, respectively, for the interval between the two tuffs. These values provide useful constraints for interpolating the fossil ages in the adjacent sediments.

The age of the Chari (= L) Tuff is less well constrained than the other tuffs owing to its large age dispersion and lower precision on individual ⁴⁰Ar/³⁹Ar analyses caused by smaller feldspar grain sizes and lower K contents for many grains (Table S5). In addition, the Chari/Tuff L samples yield a less robust Bayesian weighted mean age estimate (1359.3 ± 3.8/4.2 ka; MSWD = 3.8; P = 0.00), with the age estimate for sample 81-219 being distinctly older than those from the remaining five datasets, possibly due to the greater retention of inherited Ar. Omission of this sample from the weighted mean calculation gives a more concordant result of 1357.5 ± 1.8 ka (±2.5 ka; MSWD = 0.9; P = 0.45), which is our preferred age for the Chari (= L) Tuff.

This age is within uncertainty of previous (imprecise) K–Ar ages for the Chari Tuff (Fig. 2, Table S1), but is distinctly younger than the weighted mean age of 1386 ± 10 ka reported by McDougall and Brown (2006). This discrepancy is attributed to the different statistical approach adopted in the earlier study. Application of the Bayesian algorithm to the McDougall and Brown (2006) ⁴⁰Ar/³⁹Ar results gives a weighted mean age estimate of 1359.2 ± 7.5 ka (MSWD = 1.9), indistinguishable from the current results (Fig. 9).

The Chari (= L) Tuff and BWT in the Konso Formation have similar tuff glass chemistries (Katoh et al. 2000), but the reported BWT age of 1440 ± 40 ka is notably older than current and earlier Chari/Tuff L (McDougall and Brown 2006) determinations (Table S1). As was the case for the TRT, the youngest BWT age of 1375 ± 80 ka is within the uncertainty of the Chari/L ages, supporting correlation of the tuffs.

The more consistent Bayesian weighted mean age estimates for the three Gele Tuff samples are again preferred for this locality: 1314.9 ± 2.3 ka (±2.8 ka; MSWD = 1.9; P = 0.15) or 1315.4 ± 1.9 ka (±2.5 ka; MSWD = 0.3; P = 0.75) if the youngest age from sample TB-021-P1 is excluded (Table 1). These ages are indistinguishable from an early, imprecise K–Ar age of 1250 ± 60 ka (Brown et al. 1985), but younger than the weighted mean ⁴⁰Ar/³⁹Ar age (1331 ± 8 ka) reported by McDougall and Brown (2006);(Table 1, Fig. 2, Table S1). As before, the application of the Bayesian approach to the latter data gives a younger, analogous age (1307 ± 9 ka). The Gele and Chari tuffs do not crop out in the same sections and the current results confirm their inferred stratigraphic positions.

Uncertainties in the current tuff ages (±2–3 ka) are less than, or similar to, those typically assigned to astronomically tuned, deep-sea core sequences (±2–5 ka; e.g. Channell et al. 2020). As a result, the tuffs provide unambiguous anchor points for the astronomical tuning of adjacent sedimentary units and facilitate the testing of terrestrial palaeoclimate and palaeoenvironmental models. In turn, astronomically tuned sedimentary sequences have the potential to provide improved interpolation of fossil ages compared with estimates based on bracketing tuff and geomagnetic time markers or (often highly variable) sedimentation rates (Joordens et al. 2011, 2013).

Studies aimed at reconstructing palaeoclimate records for Turkana Basin sediments have utilized either low resolution sequence stratigraphy (e.g. Lepre et al. 2007; Nutz et al. 2017) or millennial timescale geochemical proxy records, such as hydrogen plus carbon isotopic measurements of leaf waxes (Lupien et al. 2018, 2020) and ⁸⁷Sr/⁸⁶Sr isotopic data from marine fossils (Joordens et al. 2011, 2013). Lupien et al. (2018, 2020) acquired leaf wax records from palaeolake sediments in drillcore WTK13 (see Campisano et al. 2017) located on the NW margin of Lake Turkana. The upper section of core material includes the Chari, Etirr and Ebei tuffs, as well as an unnamed tuff spanning an interval from c. 1500 to 1370 ka. However, this interval also has the lowest density of dD_wax and dC_wax isotopic data, which hindered the tuning of these proxy records (Lupien et al. 2018). By contrast, a more continuous marine fossil ⁸⁷Sr/⁸⁶Sr isotope record is available for the Upper Burgi Member (UBM), immediately below the KBS (= H2) Tuff, extending beyond the base of the Olduvai Subchron (bC2n) and possibly including a short (normal polarity) excursion thought to be the pre-Olduvai event in Area 131 (Joordens et al. 2011, 2013).

The ⁸⁷Sr/⁸⁶Sr data for fish fossils recovered from palaeolake sediments (Lake Lorenyang phase) below the KBS Tuff in the UBM of the Koobi Fora Formation (Fig. 10) have been interpreted to reflect orbitally modulated palaeoclimate (c. 21 ka precession) cycles (Joordens et al. 2011, 2013) driven by fluctuating Mediterranean monsoonal intensity (e.g. deMenocal 1995; Potts 2013). The Sr isotope palaeoclimate proxy is based on the premise (i.e. unsupported by analyses) that the palaeo-Omo River in the north drained the low ⁸⁷Sr/⁸⁶Sr (c. 0.704) volcanic rocks of the Ethiopian Highlands, whereas the palaeo-Kerio and palaeo-Turkwel rivers to the south drained metamorphic terranes characterized by higher (>0.720) ⁸⁷Sr/⁸⁶Sr isotopic values (Joordens et al. 2011, 2013). Accordingly, monsoonal maxima (wet periods) that increased drainage from the palaeo-Omo River are considered to have lowered ⁸⁷Sr/⁸⁶Sr ratios in the lake, whereas periods of reduced monsoonal strength (dry periods) led to increased lake ⁸⁷Sr/⁸⁶Sr ratios (Joordens et al. 2011, 2013; cf. van der Lubbe et al. 2017).

In an initial study, Joordens et al. (2011) tuned the UBM Sr isotopic record relative to an age of 1945 ± 4 ka (Lourens et al. 2004) for the base of the Olduvai geomagnetic Subchron (bC2n). However, astronomical age estimates for the bC2n excursion are contentious and vary by c. 43 ka (Shackleton et al. 1990; Lourens et al. 2004; Lisiecki and Ramo 2005; Channell et al. 2020), about two precession cycles, thus limiting the utility of bC2n as an anchor for orbital tuning models.

The new astronomically calibrated KBS Tuff age represents a more robust anchor point for testing the available UBM tuning models and hence the utility of the Sr isotope proxy approach (Fig. 10). A revised tuning model (Model I) for the sedimentary section below the KBS Tuff in the UBM (Areas 102 and 105; Koobi Fora Formation) was constructed using published ⁸⁷Sr/⁸⁶Sr isotopic data from marine fossils (Joordens et al. 2011, 2013), available stratigraphic logs (Joordens et al. 2011, 2013), the 65° N summer insolation model of the La2004 astronomical solution (Laskar et al. 2004) and the astronomically calibrated KBS Tuff age of 1879.1 ± 2.4 Ma (Table 1). As per previous work (Joordens et al. 2011, 2013), we assumed variable sedimentation rates across the section and initially tuned maximum Sr isotopic values to insolation minima (i.e. drier cycles) and minimum Sr isotopic values to insolation maxima (i.e. wetter cycles) (Model I) (Table S7). This tuning option is similar to Model 4 of Joordens et al. (2011). In the case of Model II, we tuned the maximum Sr isotopic values to insolation minima above the more sandy units in the UBM section (Fig. 10) relative to the above KBS age and we tuned the maximum Sr isotopic values to insolation maxima below the sandy units relative to an age of 1968 ± 2 ka for bC2n (Lisiecki and Ramo 2005) (Fig. 10; Table S7).

Our revised tuning Model I, calibrated to the KBS Tuff, includes six precession cycles that align with sapropel (Lourens et al. 1996) and dust records (Tiedemann et al. 1994) for Mediterranean and Atlantic Ocean ODP sites (Fig. 10). The alignment of the KBS Tuff on the younger side of an insolation minimum (⁸⁷Sr/⁸⁶Sr maximum) is consistent with the model assumptions (Joordens et al. 2011) (Fig. 10). However, the bC2n excursion aligns with an insolation minimum (Fig. 10) rather than the insolation maximum inferred from marine records (Lisiecki and Ramo 2005; Channell et al. 2020). In this case, the tuned ages of the bC2n and older pre-Olduvai geomagnetic transitions are 1975 ± 2 ka and >1997 ka (Fig. 9), distinctly older than previous estimates of 1968–1925 ka and c. 1977 ka, respectively (see Channell et al. 2020 and references cited therein). Furthermore, Model I contrasts with the ⁸⁷Sr/⁸⁶Sr record for Lake Turkana during Holocene times, where elevated ⁸⁷Sr/⁸⁶Sr ratios correlate with an insolation maximum (van der Lubbe et al. 2017).

These inconsistencies cast doubt on the underlying assumptions of the UBM ⁸⁷Sr/⁸⁶Sr isotope proxy tuning models. Possible reasons for this include: (1) the incorrect assignment of ⁸⁷Sr/⁸⁶Sr ratios to the palaeo-Omo v. palaeo-Kerio and palaeo-Turkwel rivers – for example, modern ⁸⁷Sr/⁸⁶Sr values for these drainages are all analogous (T. Cerling pers. comm. in van der Lubbe et al. 2017); (2) a non-continuous palaeoclimate proxy record in the UBM sediments due to the influence of local drainages and/or the reworking of sediments and fish fossils; and (3) fluctuations in local and/or regional hydrodynamic patterns related to changing palaeoclimate factors, such as the east–west migration of the Intertropical Convergence Zone and the onset of the Pacific Walker circulation at c. 2.0–2.2 Ma (van der Lubbe et al. 2017, 2021).

It is possible to reconcile these above inconsistencies by constricting the ⁸⁷Sr/⁸⁶Sr cycles by half a precession cycle (i.e. c. 10.5 ka) in the middle of the UBM succession (Model II) (Fig. 10). It is notable that this part of the stratigraphy is dominated by slightly coarser sandy units, rather than lacustrine silts (Fig. 10), and coincides with a step-change in the average ⁸⁷Sr/⁸⁶Sr ratios (Fig. 10). Although possible triggers for the shift are unclear, this period (c. 1.9–1.85 Ma) was characterized by increasing lake levels, the volcanic closure of drainage outlets to the Indian Ocean (2.2–1.8 Ma) (Boës et al. 2018 and references cited therein) and changing palaeoclimate patterns (van der Lubbe et al. 2017, 2021). Assuming that the sand units reflect a change in hydrodynamic regime, we constructed a revised tuning model (Model II, Fig. 10), with the younger record tuned to the KBS tuff age and the older record tuned to a bC2n age of 1968 ± 2 ka (Lisiecki and Ramo 2005) (Table S7). The latter age is preferred, as the younger age estimates for bC2n (1925 and 1945 ka; Horng et al. 2002; Channell et al. 2020) are incompatible with the number of definable Sr cycles and would invalidate the ⁸⁷Sr/⁸⁶Sr proxy assumptions. Although speculative, this tuning model is consistent with the ages for the KBS Tuff and bC2n (tuned age = 1965 ± 3 ka) and maintains alignment with Mediterranean sapropel (Lourens et al. 1996) and Atlantic Ocean dust (Tiedemann et al. 1994) records (Fig. 10). Nonetheless, further testing of the ⁸⁷Sr/⁸⁶Sr palaeoclimate proxy approach is clearly warranted. Our modelling demonstrates the particular challenge in tuning terrestrial palaeoclimate proxy records and highlights the importance of robust tuff (and geomagnetic excursion) ages for anchoring terrestrial tuning models in the Turkana Basin and other terrestrial settings.

Temporal constraints on hominin fossils recovered from Turkana Basin sediments rely on the ages of the bracketing tuff, geomagnetic and stratigraphic time markers, taking into account the associated uncertainties (e.g. Feibel et al. 1989, 2009; McDougall and Brown 2006). Where stratigraphic logs indicate more or less continuous deposition (i.e. no significant hiatus), hominin fossil ages have been interpolated using stratigraphic scaling (e.g. Feibel et al. 1989, 2009; Lepre and Kent 2010, 2015). Astronomically tuned sedimentary sequences, such as the UBM (Fig. 10), have the potential to provide improved age assignment for sediment/fossil deposition (Joordens et al. 2011, 2013). However, as the Sr isotope tuning model is unproven, the attribution of hominin fossil ages using this approach may be premature. The improved age precision of tuffs dated in the current study provides an opportunity to refine the temporal estimates for some associated fossils, although broader revision must await the high-precision geochronology of additional tuffs.

Relatively few hominin fossils occur in close stratigraphic proximity to the Chari (= L) or Gele tuffs. However, well-preserved early Acheulean tools occur in strata stratigraphically underlying the BWT in the Konso region (Beyene et al. 2013), which is correlated with the Chari Tuff. In addition, an unusual bone hand-axe was recovered from the same region, located c. 9 m above the BWT and assigned an age of c. 1.4–1.3 Ma (Sano et al. 2020). Our preferred age for the Chari/L/BWT Tuff revises the maximum age of the hand-axe to 1357.5 ± 2.5 ka. More significantly, the presence of the Chari Tuff in the WTK13 core offers robust anchoring of future high resolution palaeoclimate proxy data.

The c. 1.9 Ma time period that encompasses the KBS (= H2) and Malbe (= H4) tuffs is more notable archaeologically owing to the co-existence of multiple hominin species in the region, the appearance of Homo erectus, increased Homo brain sizes, the development of Acheulean stone tools and the commencement of Homo dispersion beyond Africa (e.g. Leakey et al. 2012; Potts 2013; Herries et al. 2020). This period also coincided with significant environmental changes, including elevated lake levels in eastern Africa, the expanded distribution of C4 vegetation and an upsurge in bovid fauna diversity (Levin 2015 and references cited therein). The KBS and Malbe tuffs occur in close association with Oldowan (Mode 1) stone tools (Leakey 1970), as well as several important hominin fossils, including H. erectus and early non-erectus Homo members of the ‘1813’ and ‘1470’ Groups (formerly Homo habilis and Homo rudolfensis) (Leakey 1973, 1974; Leakey et al. 2012; Anton et al. 2014).

We determined hominin fossil ages for the KBS and UBM stratigraphic sections immediately above and below the KBS and Malbe tuffs. Revised ages for hominin fossils recovered from sediments in Areas 12, 102, 105, 106, 107, 123 and 131 were determined using published stratigraphic sections (Lepre and Kent 2010; Joordens et al. 2011, 2013), a range of possible sedimentation rates, new ages for the KBS and Malbe tuffs, age estimates for the Olduvai geomagnetic transitions (tC2n and bC2n) and distinctive stratigraphic markers (mollusc-rich C4 and A2 stromatolitic markers) (Table S8). To account for the variability in tC2n and bC2n age estimates (Channell et al. 2020), median values were calculated using the full range of available estimates (1776 ± 12 ka and 1946 ± 26 ka, respectively) (Table S8). Sedimentation rates were estimated based on lithological composition and thickness, assuming essentially continuous sediment accumulation in designated sections. As sedimentation rates are typically highly variable, fossil age determinations are restricted to specimens recovered within c. 15 m of time markers (with some exceptions), where the impact of sedimentation variability is more limited. As fossils may have been reworked upwards in the sequence, these ages are essentially minimum estimates. Propagated uncertainties assigned to ages include uncertainties in the ages of the tuffs, geomagnetic and stratigraphic markers, stratigraphic position relative to chronological markers (±1 m) and the range of estimated sedimentation rates (Table S8).

The widespread C4 biostratigraphic marker occurs a few metres above the KBS and Malbe tuffs (Feibel et al. 1989) and is associated with a number of early hominin fossils in Area 123 (c. 40 km south of Area 105), including crania from the important early hominin Group 1813 (H. habilis) (e.g. Feibel et al. 2009). The timing of the C4 marker is constrained by its location c. 13 m above the KBS Tuff in section 102-17 U (Koobi Fora type section) and c. 10 m above the Malbe Tuff in section 102-18 (Feibel 1983). Based on assumed sedimentation rates of 20–30 cm ka⁻¹ for the KBS-C4 section, which contains few sandy units, we obtain an age of 1827 ± 12 ka for the C4 marker, comparable with a previous estimate of c. 1840 ka (Feibel et al. 2009). For UBM sediments above C4, which contain thicker sandy units, we assumed sedimentation rates of 20–60 cm ka⁻¹ (Table S8).

With reference to the revised C4 age of 1827 ± 12 ka, we estimated ages for Area 123 specimens located within c. 15 m of the C4 marker. Ages range from 1887 ± 47 ka (KNM-ER 1502) to 1812 ± 17 (KNM-ER1811) (Fig. 11; Table S8) and are marginally younger than prior estimates (Feibel et al. 2009), with reduced uncertainties. Revised age estimates for early Homo (Group 1470; H. rudolfensis) specimens from Area 131 (Leakey 1973; Leakey et al. 2012) range from 1984 ± 27 ka (KNM-ER 62000) to 1899 ± 7 ka (KNM-ER 1483) (Fig. 11; Table S8). Thus Group 1813 (H. habilis) specimens overlap temporally with other early Homo species, including Group 1470 (H. rudolfensis) (<1.88 to c. 2.06 Ma; Leakey et al. 2012; Joordens et al. 2013), reaffirming the presence of multiple early hominin species in the region at c. 1.9 Ma.

Previous age estimates for the important early H. erectus cranium KNM-ER 3733, the most complete example from the Koobi Fora Formation, range from c. 1780 ka (Feibel et al. 1989) to c. 1630 ka (Lepre and Kent 2015). Using the stratigraphic correlations of Tindall (1985), the Area 104 lithological section of Lepre and Kent (2015), a magnetostratigraphic age of 1776 ± 11 ka for the A2 marker bed (Area 104; Lepre and Kent 2015), a revised Bayesian age estimate for the Morte Tuff (1497 ± 6 ka; data from McDougall and Brown 2006) and a median age for tC2n (1776 ± 11 ka) (Fig. 11, Table S8), we calculate an age of 1591 ± 39 ka for KNM-ER 3733 (Table S8), marginally younger than prior estimates (McDougall et al. 2012; Lepre and Kent 2015). Older hominin fossil fragments attributed to early H. erectus, KNM-ER 2598 and KNM-ER 3228 (Leakey 1974; Wood and Leakey 2011), were recovered from normally magnetized sediments below the KBS Tuff. We obtain an age of 1915 ± 15 ka for the latter specimen, consistent with previous estimates of c. 1921 ka (Leakey 1974) and 1950 ± 50 ka (Feibel et al. 1989) (Fig. 11; Table S8). A significantly older age of c. 2040–1950 ka was reported for a H. erectus (DNH134) cranium from the Drimolen site in South Africa, based on (relatively imprecise) US-ESR molar ages of 2040 ± 48 ka (2$σ$⁠) and 1960 ± 29 ka (2$σ$⁠), a location below bC2n and the absence of the c. 2080 ka Huckleberry Ridge magnetic reversal event in the underlying strata (Herries et al. 2020). If correct, DNH134 would be the oldest known occurrence of H. erectus. However, it is notable that the c. 1977 ka pre-Olduvai geomagnetic excursion was not recognized in the Drimolen site (Herries et al. 2020) despite its identification in other Malapa cave deposits. It is therefore possible that DNH134 is located stratigraphically higher than previously thought, between bC2n and the pre-Olduvai excursion (i.e. c. 1977–1925 ka). Nonetheless, if the Drimolen Cave magnetostratigraphy is correct, then this raises the possibility that H. erectus evolved in southern Africa before dispersing northwards to eastern Africa and then Eurasia.

⁴⁰Ar/³⁹Ar dating of pumice feldspar grains from the KBS/H2, Malbe/H4, Chari/Tuff L and Gele Tuff localities in the Turkana Basin yield variably dispersed, positively skewed age distributions. The KBS/H2 and Malbe/H4 samples are characterized by relatively restricted ranges in age (36.7 and 71.7 kyr, respectively) and feldspar composition, with no relationship between the feldspar ages and composition. By contrast, the Chari/Tuff L and Gele samples display broader age ranges (>250 kyr), distinct feldspar compositional arrays and discernible correlations between feldspar ages and compositions. Consequently, the older ages are attributed to the retention of inherited Ar by feldspar antecrysts/phenocrysts in response to ‘cold storage’ in crystal reservoirs followed by rapid melt infiltration and eruption. The combination of ⁴⁰Ar/³⁹Ar age and glass/feldspar compositional data provides insights into the longevity and dynamics of (still unidentified) source magma reservoirs.

An evaluation of statistical approaches reveals that the most concordant ⁴⁰Ar/³⁹Ar eruption ages derive from the Bayesian estimation approach. This method gives preferred, astronomically calibrated ages of 1879.1 ± 0.6 ka (±2.4 ka including external errors) for the KBS (= H2) Tuff, 1837.4 ± 0.9 ka (±2.4 ka) for the Malbe (= H4) Tuff, 1357.5 ± 1.8 ka (±2.5 ka) for the Chari (= L) Tuff and 1315.4 ± 1.9 ka (±2.5 ka) for the Gele Tuff. With some exceptions, these ages are consistent with previous age determinations, but are more precise.

As the uncertainties in the current tuff ages (±2–3 ka) are similar to those of astronomically tuned deep-sea core sequences (±2–5 ka), the tuffs provide important anchor points for astronomical tuning of associated sedimentary sequences. Tuning of published ⁸⁷Sr/⁸⁶Sr proxy data from the UBM relative to the KBS Tuff age reveals some inconsistencies with geomagnetic excursions and Holocene data, suggesting complex interactions between palaeoclimate and geological drivers of sedimentation.

The new age results also permit updated age estimates for important early Homo fossils and stratigraphic marker horizons in the KBS Member and UBM. Revised ages for Group 1813 (H. habilis) and Group 1470 (H. rudolfensis) fossils range from 1887 ± 47 ka (KNM-ER 1502) to 1812 ± 17 (KNM-ER 1811) and from 1984 ± 27 ka (KNM-ER 62000) to 1899 ± 7 ka (KNM-ER 1483), respectively, consistent with multiple early hominin species occupying the Turkana Basin at c. 1.9 Ma. The calculated ages for the important H. erectus specimens KNM-ER 3733 and KNM-ER 3228 are 1591 ± 39 ka and 1915 ± 15 ka, respectively.

These results highlight the importance of a high resolution temporal framework for human evolution and palaeoclimate/palaeoenvironmental studies. Additional high-precision ⁴⁰Ar/³⁹Ar studies of Turkana Basin tuffs are needed to further test models relating palaeoclimate and palaeoenvironmental variations to hominin evolution and migration.

This study is dedicated to the memory of Ian McDougall and Frank Brown, who passed away in 2018 and 2017, respectively. This study could not have been undertaken without their encyclopaedic knowledge of the Turkana Basin. S. Szczepanski is thanked for technical assistance in the Melbourne ⁴⁰Ar/³⁹Ar laboratory. G. Hutchinson is acknowledged for assisting with the electron microprobe analyses. We also thank A. Savelkouls and H. Dalton for drafting various figures, S. Samim for assisting with sample preparation and J. Silver for generating the inverse isochron plots. The manuscript has benefitted from formal reviews by C. Hall and an anonymous reviewer.

DP: conceptualization (lead), data curation (lead), formal analysis (lead), funding acquisition (lead), investigation (lead), methodology (lead), project administration (lead), resources (lead), supervision (lead), validation (lead), visualization (lead), writing – original draft (lead), writing – review and editing (lead); ELM: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), methodology (equal), project administration (equal), resources (supporting), supervision (supporting), validation (supporting), writing – original draft (equal), writing – review and editing (equal); AJG: conceptualization (supporting), funding acquisition (supporting), investigation (supporting), resources (supporting), supervision (supporting), validation (supporting), writing – original draft (supporting), writing – review and editing (supporting); FB: conceptualization (supporting), investigation (supporting), methodology (supporting), project administration (supporting), resources (supporting), writing – original draft (supporting); IM: conceptualization (supporting), investigation (supporting), methodology (supporting), resources (supporting), supervision (supporting), writing – original draft (supporting); TEC: methodology (supporting), resources (supporting), supervision (supporting), validation (supporting), writing – original draft (supporting), writing – review and editing (supporting); ML: conceptualization (supporting), funding acquisition (supporting), resources (supporting), supervision (supporting), writing – original draft (supporting), writing – review and editing (supporting); JH: formal analysis (supporting), funding acquisition (supporting), methodology (supporting), resources (supporting), supervision (supporting), writing – original draft (supporting), writing – review and editing (supporting); LL: conceptualization (supporting), project administration (supporting), resources (supporting), validation (supporting), writing – original draft (supporting), writing – review and editing (supporting).

This study was supported financially by Australian Research Council Discovery Grant DP180101412 to DP, AJWG, JH, ELM, IM and ML. The Melbourne ⁴⁰Ar/³⁹Ar laboratory receives ongoing operational support under the AuScope programme of the National Collaborative Research Infrastructure Strategy (https://www.education.gov.au/ncris; https://www.auscope.org.au).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. To the best of their ability, the Journal has determined that Professors Brown and McDougall worked on this project and that there were no competing interests.

All data that support this study research are included in the article and/or supporting information.