Integrated astrochronology, sequence stratigraphy, and chronostratigraphy of a shallow marine sandy mudstone (Lower Jurassic, Redcar Mudstone Formation, Cleveland Basin, UK)

Abstract

Cyclostratigraphy and sequence stratigraphy are cornerstone disciplines of modern sedimentary geology but are not commonly applied to the same study sections. Here, we apply both concepts to the argillaceous, shallow marine, Early Jurassic (Sinemurian–Pliensbachian) Redcar Mudstone Formation exposed at Robin Hood’s Bay, Yorkshire. Hand-held X-Ray Fluorescence (HH-XRF) high-resolution elemental data and bulk organic matter carbon-isotope (δ¹³C_org) data are presented to elucidate the context and pace of climatic change in the Early Jurassic, focussing on the Sinemurian–Pliensbachian boundary interval. The carbon-isotope stratigraphy demonstrates a 4.5 ‰ negative carbon-isotope excursion in bulk organic matter, representing the Sinemurian–Pliensbachian Boundary Event (SPBE). Elemental ratios (Si/Al and Zr/Rb), which are effective grain-size proxies, reveal coarsening upward cycles at a scale of 1.2–1.6 m, and are interpreted as parasequences. Parasequence stacking patterns allow the identification of depositional sequences, signifying five cycles of relative sea level change, with a major 2nd order relative sea-level rise and transgression associated with the Sinemurian–Pliensbachian boundary event. Missing or condensed retrogradational parasequences are associated with flooding and indicative of hiatuses at sub-biostratigraphic resolution – particularly in the Oxynotum Subzone of the Oxynotum Zone, and the Macdonnelli Subzone of the Raricostatum Zone. Maxima of redox-sensitive elements (Mo, V) and TOC in the Jamesoni Zone indicate dysoxic palaeoredox conditions in an otherwise oxic succession at the start of the Pliensbachian. Spectral analysis indicates regular cyclicity for parasequences through the succession, which are interpreted as being forced by short orbital eccentricity, and abrupt transitions in evolutive harmonic analysis support the hypothesis of short-duration hiatuses identified using sequence stratigraphy. The estimation of missing time is constrained using the differences in durations from models tuned to interpreted 100 kyr and 405 kyr orbital eccentricity cycles. The tuned age model gives estimates of the duration of the Oxynotum (0.90–1.00 Myr) and Raricostatum (1.35–1.45 Myr) zones. This further suggests a duration of 1.1 Myr for the early phase of the Sinemurian–Pliensbachian boundary event and 200 kyr for the Sinemurian age Liasidium Event, both regarded as possible hyperthermal episodes. Integrated sequence stratigraphy and astrochronology help to refine the timescale for the Early Jurassic and provide a framework for understanding the precise timing of the development of sequence stratigraphic surfaces.

Lay summary

The amount of sunlight received at the Earth’s surface is strongly affected by periodic cycles in planetary orbits and, in turn, these cyclical changes are recorded in the sedimentary rock record. Identification of cycles can thus be used to provide estimates for the time taken to deposit sedimentary strata. However, the same climatic parameters also affect sea-level changes, and therefore, completeness of the sedimentary record. Here, we combine the study of both small-scale lithological climatic cycles (Milankovitch Cycles) and larger-scale cycles forced by sea-level change (Sequence Stratigraphy) in order to estimate durations for Early Jurassic time increments (Chronozones) at an international reference location in the Cleveland Sedimentary Basin, Yorkshire, UK.

1. Introduction

Cyclostratigraphy and sequence stratigraphy are powerful tools to understand the sedimentary record (Laskar 2020; Simmons et al., 2020; and references therein). Cyclostratigraphy refers to the astronomical forcing of climate and its expression as quasi-periodic cycles within sedimentary successions (see Weedon, 2022, for application to Jurassic stratigraphy). Cyclostratigraphy, in conjunction with biostratigraphy, magnetostratigraphy and radioisotopic dating, has the power to establish age models tuned to Earth’s astronomical parameters (e.g., Hays et al., 1976; Imbrie et al., 1984; Laskar, 2020; Weedon, 2003), and is thus also known as astrochronology.

Sequence stratigraphy was developed as a framework to interpret the temporal-spatial distribution of shallow marine sediments in response to palaeoenvironmental and geological processes (Vail, 1977; Posamentier & Vail, 1988; Posamentier et al., 1988). Sequence stratigraphy is built on a hierarchical system of sedimentary units: parasequences, systems tracts and sequences (Posamentier et al., 1988; Van Wagoner et al., 1988). Each unit is bounded by time-significant surfaces related to non-deposition or erosion (Sloss, 1963; Posamentier et al., 1988; Van Wagoner et al., 1988). In shallow marine siliciclastic settings, sequence stratigraphic surfaces and internal facies arrangements, such as coarsening and fining upwards trends, are strongly influenced by relative sea-level change; a term that encompasses the mixed effects on basin-scale sedimentary accumulation of sediment loading, tectonic activity and eustatic sea level (Milne et al., 2009). The full spectrum of processes that influence the sequence stratigraphic record was recently summarised and reviewed by Bohacs (2022).

There clearly must be a relationship between astronomical forcing and sequence stratigraphy because two of the principal factors that determine sequence architecture (eustatic sea-level change and variation in sediment supply) are also affected by the astronomically controlled variations in climate (Miller et al., 2020; Simmons et al., 2020). Conversely, the predictable geometries of depositional sequences have implications for sedimentation rate and stratigraphic completeness, knowledge of which is essential for the construction of robust age models and timescales. Furthermore, the integration of cyclostratigraphy and sequence stratigraphy offers a route to understand the extent to which sequence stratigraphic expression in the rock record is determined by eustatic sea-level change, which requires synchronous development across multiple basins (e.g., Haq, 2017) or other more localised and complex tectonic or sedimentary processes (e.g., Willis et al., 2022).

1.1. Background and aims

The Early Jurassic shows evidence of two major hyperthermal events, at the start and near the end of the epoch, namely the Triassic-Jurassic boundary event 201 Myr ago (e.g., Pálfy et al., 2001; Korte et al., 2019) and early Toarcian Oceanic Anoxic Event 183 Myr ago (e.g., Jenkyns, 1985, 1988; Thibault et al., 2018; Ruhl et al., 2022). Research into the ~18 million years between these big events has drawn attention to several smaller magnitude perturbations to the carbon cycle, climate and environment, which have no certain causes (Korte & Hesselbo, 2011; Riding et al., 2013; Franceschi et al., 2014; Price et al., 2016; Hesselbo et al., 2020a). In particular, it is suggested that other hyperthermal events also occurred, during the Oxynotum Ammonite Chronozone in the Sinemurian (called the Liasidium Event after a distinctive dinoflagellate taxon common in this interval in NW Europe), and at the Sinemurian–Pliensbachian boundary (Raricostatum–Jamesoni zones; Korte & Hesselbo, 2011; Riding et al., 2013; Hesselbo et al., 2020a; Franceschi et al., 2022). Note, Jurassic time is conventionally subdivided using ammonite-based chronozones (Cox, 1990; Page, 2017; Hesselbo et al., 2020b).

These two possible hyperthermals are understudied, but they both show 2–4.5 ‰ negative carbon-isotope excursions (CIE’s) recorded in bulk organic matter and bulk carbonate, with the Sinemurian–Pliensbachian boundary event also showing parallel excursions in fossil wood, bulk carbonate, and macrofossil carbonate (van de Schootbrugge et al., 2005; Suan et al., 2010; Korte & Hesselbo, 2011; Riding et al., 2013; Peti et al., 2017; Ullmann et al., 2021). Knowledge of the causes and durations of these events is restricted by a sparsity of detailed studies, as well as a lack of robust astronomically-tuned age models. Additionally, many Lower Jurassic shallow marine successions are prone to stratigraphic gaps (i.e., times of erosion and non-deposition) with few opportunities to create continuous astronomical timescales.

Milankovitch frequency cyclicity has been observed in studies of sections of Sinemurian mudstones (van Buchem & McCave, 1989; van Buchem et al., 1992; Ruhl et al., 2016; Xu et al., 2016; Storm et al., 2020). Despite this, the Sinemurian stage is currently one of the few intervals of time in the Mesozoic with no complete and resolved astronomical timescale (Leu et al., 2024), while considerable debate still surrounds the floating timescales for the Hettangian (Hüsing et al., 2014; Ruhl et al., 2016; Weedon et al., 2019) and Toarcian stages (Suan et al., 2008; Boulila et al., 2014; Huang, 2018; Hesselbo et al., 2020b).

In the context of fine-grained siliciclastic successions, an approach using relatively large geochemical datasets obtained by core scanning or handheld X-Ray Fluorescence (XRF) has recently been shown to be effective for sequence stratigraphic interpretations (e.g., LaGrange et al., 2020; Thöle et al., 2020), and as a basis for cycle analysis (e.g., Ruhl et al., 2016; Thibault et al., 2018; Hollaar et al., 2023). This paper aims to reconstruct detailed sequence stratigraphy from field observations and elemental and isotope geochemical time series in largely argillaceous, shallow-marine sedimentary strata, and integrate this approach with astrochronology to determine the pace and forcing of climate change in the Jurassic. These aims are applied to the internationally important Sinemurian–Pliensbachian boundary Global Stratotype Section and Point (GSSP) site at Robin Hood’s Bay, Yorkshire, UK (Figure 1; Hesselbo et al., 2000b; Meister et al., 2003, 2016), enabling us to anchor our new interpretation to the international stratigraphic framework.

1.2. Geological setting

Robin Hood’s Bay is located in the Cleveland Basin in North Yorkshire. The Cleveland Basin is a relatively small Mesozoic basin bounded to the north-east, west and south by the Mid North Sea, Pennine and Market Weighton structural highs (Hemingway, 1974; Powell, 2010). The Lower Jurassic comprises ~450 m of shallow marine, mostly siliciclastic deposits ranging from Hettangian to Toarcian in age (Hesselbo & Jenkyns 1995; Hesselbo & King, 2019; Simms et al., 2004). The Hettangian to lower Pliensbachian is represented by the Redcar Mudstone Formation (~300 m thick). This formation comprises beds of bioturbated, fossiliferous silty mudstone and well-cemented muddy sandstone, deposited in environments that range from the offshore transition to shoreface under storm-dominated conditions (van Buchem et al., 1992; Powell, 2010). Broadly, the upper Sinemurian is represented by the Siliceous Shale Member, with the lower Pliensbachian represented by the Pyritous Shale Member and the Ironstone Shale Member (Figure 1; Hesselbo and Jenkyns, 1995; Howarth, 2002; Hesselbo & King 2019). The Siliceous Shale is coarser-grained than the Pyritous Shale or the Ironstone Shale, and all these members have well-developed marker beds of decimetre-thick diagenetic siderite nodules (Figure 2).

The Cleveland Basin formed part of a network of extensional Mesozoic depocentres that developed across much of present-day, north-west Europe and the adjoining North American and East Greenland margins. These Early Jurassic basins linked up as an extensive shallow-marine realm known collectively as the Laurasian Seaway (Ziegler, 1990; Bjerrum et al., 2001) sometimes also called the European Epicontinental Shelf (e.g., Ruvalcaba Baroni 2018). The Cleveland Basin is thought to have received sediment from Paleozoic basement with sediment transport from north to south (Hemingway, 1974; Hesselbo and Jenkyns, 1995; Powell, 2010, and references therein).

↗

Figure 1. (A) Summary map of Robin Hood’s Bay adapted from Howarth (2002), showing the foreshore and cliff exposures. Red lines mark the positions of sampling profiles. Foreshore ammonite subchonozone abbreviations: Sau = Sauzeanum, Bro = Brooki, Bir = Birchi, Obt = Obtusum, Ste = Stellare , Dnt = Denotatus , Sim = Simpsoni, Oxy = Oxynotum, Den = Densinodulum, Rar = Raricosta toides , Mcd = Macdonelli, Apl = Aplanatum. (B) Global palaeogeography for the Early Jurassic, modified from Dera et al. (2009). The Cleveland Basin sits in an epeiric sea setting within the N-S oriented Laurasian Seaway. (C) Field photograph of Wine Haven, Robin Hoods Bay. Looking SW. Figures for scale. (D) Lithostratigraphy and biostratigraphy for the Early Jurassic at Robin Hood’s Bay (Hesselbo & Jenkyns, 1995; Howarth, 2002). Timescale from Hesselbo et al. (2020b). Ammonite zones and subzones shown are only those relevant to this study. (E) The Cleveland Basin and other referenced localities in their local geographic context (adapted from Cox et al., 1999). (F) Key to lithological log shown in Figure 2. Note that Jurassic ammonite chronozones are by convention labelled according to a species name in roman font with an initial capital letter (see Page, 2017, Hesselbo et al.,2020b, and references therein).

2. Methodology

2.1. Field data

New samples were collected at Robin Hood’s Bay and correlated to previously compiled measured outcrops based on centimetre-by-centimetre visual description (Figure 2; Hesselbo and Jenkyns, 1995; Howarth, 2002; Hesselbo et al., 2020a). Heights are recorded with respect to a level near the top of the Siliceous Shale Member, a cemented sandstone that forms a protruding ledge known locally as East Scar.

Samples were collected along a composite section, primarily from the base of the cliff, utilising the broad anticline to access progressively higher strata on the south side of the bay (Figures 1, 2). Sampling of the upper Sinemurian and Pliensbachian (between -7 m and +16 m height) was conducted at Wine Haven between [54.40652°N, 0.49350°W] and [54.40678°N, 0.49862°W], shown in Figure 1, during two periods of fieldwork in 2017 and 2018. A vertical resolution of better than 5 cm was achieved for this sample set of 500 samples. This sample set was combined with 205 samples previously collected at 10 cm vertical resolution (Hesselbo et al., 2020a), from the middle Sinemurian Obtusum and Oxynotum zones in the cliffs and foreshore between Stoupe Beck and Boggle Hole in Robin Hood’s Bay (between -7 and -28 m; Hesselbo et al., 2020a). All samples were chosen to be representative of the background lithology, and diagenetic features (such as siderite nodules) and prominent fossils were avoided unless they formed major components of continuous beds.

↗

Figure 2. Lithological log adapted and extended from Hesselbo and Jenkyns (1995) and Hesselbo et al. (2020a) ranging through the Obtusum–Jamesoni zones. Bed numbers Ref. A are from Howarth (2002); Bed numbers Ref. B are from Hesselbo and Jenkyns (1995). See Figure 1 for graphic log key. Ammonite determinations and ranges largely follow the records of Bairstow as revised by Howarth (2002). Where changes have been made in the present study (e.g., re-identification , extension of range, etc.) these are indicated by “*” and explained in the main text. Closed circles indicate occurrences reported to precise stratigraphic levels, whilst open circles indicate occurrences reported from the equivalent bed of Howarth (2002; i.e., Ref. A). Newly reported ammonite occurrences are indicated with square symbols. Ammonite horizons and biohorizons are described and summarised in Page (2003, 2010, 2017). A vector graphics version of this file with individual sample horizons shown to 1 cm precision is available as Supplementary Data File 1.

2.2. Hand-held X-Ray Fluorescence analysis (HH-XRF)

Samples were cleaned with de-ionised water and were then thoroughly dried in an oven overnight at 40°C (to avoid any moisture-related attenuation of X-rays during XRF analysis; Gazley & Fisher, 2014; Schneider et al., 2016). All samples were ground using an agate pestle and mortar to a fine homogenous powder and then placed in airtight plastic tubes.

HH-XRF analysis was conducted on powdered bulk rock samples using two Olympus Innov-X Delta Premium XRF analysers at the University of Exeter Camborne School of Mines. The devices were operated in Geochem mode and readings were collected for 120 s each at voltages of 40 kV and 10 kV (240 s total analysis time). Elements selected for analysis include Si, Al, Zr, Rb, Mo, V, Ni and S. Samples were placed in 30 ml plastic centrifuge tubes to obtain a depth of at least 1 cm of compacted sediment. The tubes were covered with low density polyethylene film of ~10 µm thickness, inverted and placed into the analysis chamber of the Olympus XRF stand. Reproducibility was controlled through a total of 56 interspersed measurements of the standard NIST 2711a (dried, powdered soil from Montana, USA, with moderately elevated trace element concentrations). Due to matrix and instrumental effects, HH-XRF results can only be regarded as semi-quantitative.

Nevertheless, the relative trends in the readings are thought to be robust as evidenced by standard reproducibility. Reproducibility of the NIST standard (2 relative standard deviations) is 6 % for Si/Al, 4 % for Zr/Rb, 13 % for S/Fe, 15 % for V, and 26 % for Ni. Detection limits of the technique with the chosen settings are irrelevant for Al, Si, Zr, Rb, Fe, V, and Ni, for which all measured values are well above this limit, but most samples show Mo concentrations below the detection limit. Mo concentrations in NIST 2711a are consistently below the detection limit of the instrument using the chosen settings, so reproducibility for this element cannot be quantified in this way. Instrumental estimates of detection limits (LOD) are unreliable and often several hundred % higher than the lowest quantified values. For Mo, quantification limits are likely better than 10 µg/g, but intervals where only a few samples are indicated at Mo levels of 1 to 10 µg/g should be taken with caution.

The samples were analysed using two different devices due to unscheduled maintenance on device A. Device B was the same model and specification; however, a device offset was encountered that required minor correction. A correction equation was created for each element of interest using repeat analysis of a series of representative samples from the set analysed on both devices. Linear trends arising from cross-matching the analysed values were used to correct the samples analysed on device A to match those of device B. Correlation coefficients (r²) for these linear fits were usually high, >0.99 for Rb and Fe, >0.9 for Zr and Ni, 0.82 for Si and lower for elements which showed only a moderate spread of concentrations in the calibration dataset (0.57 for Al, 0.53 for V, and 0.3 for S). For full treatment, see Supplementary Data File 2.

2.3. Bulk organic carbon isotopes and total organic carbon

Organic carbon isotope ratio (δ¹³C_org) and total organic carbon (TOC) contents were determined at the University of Exeter Environment and Sustainability Institute (ESI) stable isotope facility on the Penryn Campus, Cornwall. The analysis was conducted with a coupled Sercon elemental analyser - continuous flow - isotope ratio mass spectrometer (EA-CF-IRMS). Splits were used for samples prepared for HH-XRF analysis. The organic carbon fraction of the powdered rock sample was isolated via de-carbonation to remove any inorganic carbon associated with carbonate minerals. De-carbonation was achieved with an acid digestion method using 3.3 % hydrochloric acid heated to 70°C for 1 hour and left to cool for 24 hours. Reacted remains were washed repeatedly with de-ionised water to remove chlorides before further processing for analysis. A total weight of decarbonated sediment containing ~100 µg of organic matter was measured, placed into sealed tin capsules and loaded into the CF-IRMS auto-sampler. In-house standards Sigma Alanine and Bovine Liver that were previously calibrated to certified reference materials IAEA-N1, IAEA-N2, L-SVEC, IAEA-CH6 and B2155 were used to correct for instrument drift and to determine the δ¹³C values of the samples (Supplementary Data File 2). Total organic carbon was determined using the CO₂ beam area relative to the in-house standard. The 2 * standard deviation for the samples analysed was 0.09 ‰ for δ¹³C, and uncertainty of TOC analyses expressed as 2 * relative standard devitation is 1.1 %, based on in-house standard measurements of 47.50 +/- 0.53%.

2.4. Power spectral analysis

Analysis was conducted primarily in Acycle (v2.1; Li et al., 2019). Power spectral analysis was conducted with the 2π multi-taper method (MTM; Thomson et al., 1982) with AR-1 (auto-regressive) robust red noise modelling (Mann & Lees, 1996) to estimate statistical confidence relative to the null hypothesis of no astronomical forcing. Confidence limits are drawn at 99.9%, 99%, 95% and 90% confidence, and peaks are taken as significant if exceeding the 95% level. Detrending was conducted primarily using the Locally Estimated Scatterplot Smoothing (LOESS) method in the Acycle software package to remove unwanted long-term stratigraphic trends. Tuning to construct a floating astronomical timescale was conducted on interpreted 405 kyr and 100 kyr stratigraphic cycles through extraction using elliptic or Gaussian filters provided in Acycle.

3. High-resolution ammonite-based chronostratigraphy

For Late Sinemurian biochronology, Robin Hood’s Bay provides one of the most complete, permanently exposed successions known, and is thus fundamental to the interpretation of a significant number of ammonite species of stratigraphical importance. The location also provides actual or potential stratotypes, which facilitate the establishment of true chronostratigraphical meaning for correlative units derived from the ranges of these ammonite species. Significant biostratigraphical records are scattered through many works from the later 19^th century into the 21^st century (e.g., Tate & Blake, 1876; Wright, 1878– 1882; Buckman, 1909–1930; Spath, 1925–1926; Getty, 1973; Page, 1992; Dommergues et al., 1994; Simms et al., 2004; Meister et al., 2006).

A notable synthesis is provided by Howarth (2002), which integrates the historical sections and collections of L. Bairstow and places them into the lithological succession. Here, we update Howarth’s account by recognizing significant results from elsewhere, including at a high resolution, sub-subzonal level. The high-resolution biohorizon scheme applied here follows that developed by Corna (1985), Dommergues and Meister (1989), Page (1992), and Dommergues et al. (1994) as reviewed by Corna et al. (1997) and summarised in Page (2003). See Page (2017) for a discussion of the zonal terminology used, including biohorizons, horizons and chronozones. The new framework (Figure 2), is based on a re-assessment of the determinations and stratigraphical conclusions of Howarth (2002), combined with new observations by co-author KNP (1980s–1990s). Further refinement may in the future be possible based on Bairstow’s original specimens in the Natural History Museum, London. The principal observations and arguments that underpin the revised biochronology are summarised in the subsections below.

3.1. Obtusum Chronozone

Stellare Subchronozone — the record of Epophoceras in beds 446.4 (1 m above base) in Howarth (2002) indicates biohorizon Sn 36 (Index: ‘Galaticeras’ sp. nov.) at the base of the Stellare Subchronozone (i.e., lower than placed by Howarth), although the immediately succeeding Sn37–Sn40 are not currently well characterised. The range of Asteroceras grp blakei Spath without Aegasteroceras from beds 452 to 454.1, however, indicates Sn41 (blakei) to Sn42 (aff. arnouldi (Dumortier)), and Aegasteroceras grp sagittarium (Blake) without Eparietites from beds 454.1 to 456 indicates Sn43 (sagittarium). Note that ‘Eparietites’ bairstowi in Howarth (2002) from the upper part of this interval is interpreted here as a macroconch of a late member of the As. blakei group.

Denotatus Subchronozone — biohorizon Sn44 (cf. undaries) with Eparietites cf. undaries (Quenstedt) from beds 457 to 458.2 marks the base of subchronozone (after Page, 1992), slightly higher than placed by Howarth (2002), with Eparietites, including Ep. denotatus (Simpson) (KNP pers. obs.), from beds 463 to 464.32, Sn46 (denotatus) to Sn47 (glaber). Note that this interval was placed in the Simpsoni Subchronozone by Howarth (2002), based on the simpsoni-like smooth body chambers that macroconchs of late species of Eparietites possess (as noted by Page, 1992). The parallel record of ribbed ‘Epareititesimpendens (Simpson)’ (i.e., inner whorls/ microconchs) and ‘Oxynoticerassimpsoni’ (i.e., smooth, mature macroconch body chambers) in Howarth’s section proves Denotatus Subchronozone. However, the biohorizon Sn45 (fowleri) is not proven in Robin Hood’s Bay.

3.2. Oxynotum Chronozone

Simpsoni Subchronozone —the base of the subchronozone is here revised to the first level recorded by Howarth where ‘Ox. simpsoni (Simpson)’ occurs without associated ‘Ep. impendens’ in bed 464.33. The range of Gagaticeras sensu stricto (e.g., grps. exortum (Simpson) and gagateum (Young and Bird) from beds 467 to 470 indicates the Sn48 (exortum) to Sn49 (gagateum) biohorizons. The succeeding driani Biohorizon (Sn50) is not recognised.

Oxynotum Subchronozone — the base of the subchronozone can be correlated with the first occurrence of Ox. grp oxynotum (Quenstedt) in bed 472.1, which indicates Sn51 (grp. oxynotum). Bifericeras sp. in beds 482.4 (KNP obs.) to 483.1 indicates Sn52 (bifer), and the interval above the last recorded Oxynoticeras with Gleviceras only (beds 484 to 485.2) indicates Sn53 (doris).

3.3. Raricostatum Chronozone

Densinodulum Subchronozone — the basal Biohorizon of the subchronozone, Sn54 (delicatum) and the corresponding Delicatum Horizon are not confirmed in the Bay, but if present, would be expected for beds 485.2 (i.e., immediately above Glevicerasdoris (Reynès)) and 486.1, because ‘Microceras’ grp subplanicosta (Oppel) is present in 486.2 (recorded by Howarth [2002] as ‘Bifericeras cf. vitreum (Simpson)’) indicating Sn55 (subplanicosta). Crucilobiceras densinodulum Buckman in beds 486.3–487 indicates a level very close to Sn56 (lymense) and echioceratids in beds 488–489 potentially signify Sn59 (radiatum) and Sn60 (‘Echioceras sp. 3’) at the top of the subchronozone, although the intervening Sn57 (bispinigerum) and Sn58 (grp. ‘armatum’) have not been detected.

Raricostatoides Subchronozone —the presence of raricostate echioceratids from near the base of 491 (KNP obs.) to 491.1 suggests the early part of the subchronozone, including Sn61 (rhodanicum), Sn62 (raricostatum) and Sn63 (crassicostatum) with Echioceras grp intermedium (Trueman and Williams) from beds 491.2 to 493.2 indicating Sn64 (cf. intermedium) and Sn65 (boehmi).

Macdonnelli Subchronozone —following Page (1992) and Dommergues et al. (1994) the basal biohorizon of the subchronzone is taken as Sn66 (subplicatum), as indicated by Leptechioceras planum Trueman and Williams from beds 493.3 to 493.5. This level was formerly included in the Raricostatoides Subchronozone by Howarth (2002).

Aplanatum Subchronozone — Howarth’s (2002) record of ‘Paltechiocerasregustatum (Buckman)’ in 496 indicates the base of the subchronozone and probably Sn69 (aureolum), which, according to the former’s interpretation, ranges up to a level in bed 497. ‘Pa. tardecrescens (Hauer)’, similarly from bed 497 up to the basal c. 10 cm of 499, indicates both Sn70 (recticostatum) and Sn71 (tardecrescens).

3.4. Jamesoni Chronozone

The first occurrence of ‘Bifericeras’ donovani Domergues and Meister in Bed 501.1 corresponds to the defined base of the Pliensbachian Stage (Meister et al., 2006) and the first subzone of the Jamesoni Zone (Taylori Subzone). When the GSSP was ratified, this was a unique record of the indicator species (at its type locality), but it is now supplanted by records from both the Llanbedr (Mochras Farm) borehole in NW Wales, and the Prees 2 borehole in Shropshire (Hesselbo et al., 2023), hence confirming its utility as a basal Pliensbachian indicator. At Robin Hood’s Bay, as is the case for the Mochras, Prees and Stowell Park (Gloucestershire; Spath, 1956) cores, there is a Sinemurian-Pliensbachian ‘Tubellites’ tubellus (Simpson) fauna (assigned by Howarth [2002] to ‘Gemellaroceras’) that spans the boundary from at least beds 497 to 505.1.

4. Results

4.1. Organic matter and organic carbon isotopes

An extended organic carbon-isotope stratigraphy from bulk organic matter for the Cleveland Basin is shown in Figure 3. The new 10 cm resolution carbon-isotope data series extends upwards from the Sinemurian Oxynotum Zone to the mid Taylori Subzone (Jamesoni Zone) of the Pliensbachian Stage to include the base-Pliensbachian GSSP, complementing previously published data from Hesselbo et al. (2020a), which is also at 10 cm resolution for the lower part of the succession. Note that HH-XRF data were generated at full 5 cm resolution for the upper sample set.

↗

Figure 3. δ¹³C_org values and Total Organic Carbon (TOC) contents for the Robin Hood’s Bay section. Visually determined bundles of cycles (hypothesised to represent the 405 kyr long orbital eccentricity) are marked with grey dashed arcs. Intervals of elevated TOC occur from 9 m height to the top of the studied section, and negative excursions in δ¹³C_org occur from -21 m to -17 m and 9 m to 15 m height. Simplified graphic log and marker beds are based on Figure 2. Redox-sensitive elements and element ratios are also shown. Solid horizontal lines represent ammonite chronozone boundaries. Dashed horizontal lines mark the ammonite subchronozone boundaries. Hgt = height. For full ammonite chronozone names see Figure 2.

The dataset shows a significant negative carbon-isotope trend that begins in the mid Raricostatum Zone at -3 m height with a gradual (over 12 m) decrease in δ¹³C_org from -24 ‰ to -25.5 ‰ to the Sinemurian-Pliensbachian (S-P) boundary (+9 m height). The overall gradual decrease is stepped with a 6- to 4-m-thick motif of isotopically lighter values correlating with coarser-grained units as seen in the graphic log. At the Sinemurian–Pliensbachian boundary, a sharp drop is recorded from -25.5 ‰ to ~-28 ‰ (Figure 3). This negative CIE peak itself comprises roughly two 1.5-m-thick, 1 ‰ ‘spikey’ excursions, where δ¹³C_org jumps and fluctuates between -29 ‰ and -27.5 ‰ for several metres. The values do not fully recover in the measured dataset, but they increase again from ~13 m height to the top of the studied section at 16 m. Apart from a localised spike associated with the Double Band bed at ~-15 m height, the least negative values occur at bed 59 (-3.35 m), in the lower Raricostatoides Subzone.

The TOC is relatively low and stable for much of the section, particularly below 9 m height, and the lowest TOC values correspond to the relatively coarser-grained beds (Figure 3). The data show a small range between 0.2 wt% and 1.3 wt%. At 10 m height, in the lowest Pliensbachian, the TOC increases sharply from 0.75 wt% to >3 wt%. This relatively organic-rich interval is associated with the isotopically lighter carbon isotope values. Previously published organic matter pyrolysis data from Robin Hood’s Bay (van Buchem et al., 1992) show that Hydrogen Index jumps from generally less than 50 in the Siliceous Shale to around 200 in the Pyritous Shale; thus, there is undoubtedly a higher proportion of marine organic matter included in the analyses for the Pyritous Shale.

4.2. Elemental geochemistry

4.2.1. Sediment composition and grain size proxies and parasequences

For the most part, this section is siliciclastic with negligible carbonate in the rock matrix (Figure 4; see Hesselbo et al., 2020a, for petrographic characterisation of the Siliceous Shale). Calcium is mostly below detection limit but fluctuates around 0.5% from -16 to -28 m, and 0.5 to 2% from 9 to -16 m. Increases in calcium occur abruptly in association with the hard, carbonate-cemented siltstone and sandstone beds (<4%). There are also significant concentrations of carbonate-shelled fossils and nodules on the tops of sandstone beds.

The ratio of silicon/aluminium (Si/Al) varies in accordance with the grain size shown in the lithological log based on centimetre-by-centimetre visual observation in the field (Figure 4) and is also concordant with changing ichnofacies characteristic of coarser- and finer-grained depositional substrates (Sellwood, 1970; McIlroy, 2007; Powell, 2010). High Si/Al ratios correlate well with coarser-grained siliciclastic beds and lower Si/Al ratios correlate with finer-grained mudstone beds (Figure 4). This association suggests that the Si/Al ratio acts as a reliable grain-size proxy in the Robin Hood’s Bay section. Silicon is likely tracking the abundance of coarser quartz grains, an interpretation that is supported by spot-sample X-Ray Diffraction data and SEM analysis from -7 m to -28 m in the section (Hesselbo et al., 2020a). Aluminium is sourced primarily from the finer aluminosilicate (i.e. clay) minerals (Calvert & Pederson, 2007). The zirconium/rubidium (Zr/Rb) grain-size proxy shows nearly identical stratigraphic trends to the Si/Al proxy, indicating that this is also following the average grain size of the sediments. Zirconium is almost exclusively found in silt to fine sand-grade heavy sedimentary zircons and, due to its density, is associated with fine to medium sand-grade sediments (e.g., Dypvik & Harris, 2001; Calvert & Pedersen, 2007; Playter et al., 2018). Rubidium is substituted for potassium in aluminosilicates and feldspars (Heier & Adams, 1964; Playter et al., 2018). Association with fine-grained units, and generally low feldspar contents, shown in previous work, suggests Rb is mostly present in clays (Hesselbo et al., 2020a).

The utility of using Zr/Rb and Si/Al ratios as proxies for grain size has now been well established for similar depositional settings, with recent stratigraphically comparable examples found in Thibault et al. (2018), Xu et al. (2018), Thöle et al. (2020) and Hesselbo et al. (2023). Of particular relevance for the present study is the work of van Buchem et al. (1994), who analysed in detail the relationship between Si/Al ratio, quartz/clay ratio, magnetic susceptibility, and grain size for the Pliensbachian Ironstone Shale Member of the Lower Jurassic Redcar Mudstone Formation in Robin Hood’s Bay. Although in the Ironstone Shale the grain size variations are more subtle than is the case for the Siliceous Shale and Pyritous Shale members that are the focus here, a very close correspondence was found between Si/Al ratio and grain size (van Buchem et al., 1994). Of course, a diagenetic siliceous component may also potentially be present; however, although there is some limited silicification for some British Lower Jurassic successions (e.g., South Wales, Hebrides, see Korte et al., 2015; Hesselbo et al., 2020c), silicification has not been reported for the Siliceous Shale/Pyritous Shale interval. Similarly, there are no radiolarians reported from any British Lower Jurassic rocks either onshore or offshore, and no reports of other siliceous fossils from the Lower Jurassic in Yorkshire (see Lord, 2019). Furthermore, analysis of Si/Al ratios together with Zr/Rb ratios, would be effective to discover these phenomena if they were important, as also discussed by Thöle et al. (2020).

Against depth, grain size proxies show two large-scale cycles in the data. From the silty Lower Triplet bed in the Obtusum Zone, the sequence is fining upwards, culminating in the most argillaceous facies of the Sinemurian within the Simpsoni Subzone of the Oxynotum Zone. Overlying this, in the Oxynotum Subzone, several beds show coarsening upward trends that culminate in the coarsest bed, the Double Band. On the largest scale, the trend from here to the Taylori Subzone and Jamesoni Zone in the Pliensbachian (-15.5 m to 10 m), is one of gradual fining based on the Si/Al ratio, whereas Zr/Rb presents some fluctuation around a stable mean. Whether or not the overall fining begins at the Double Band or the Landing Scar level (Figure 4), this decrease in grain size reaches its minimum in the Pyritous Shale Member of the basal Pliensbachian.

Imposed on this overall trend, the grain size proxies show a saw-tooth variation on the order of 1.2 m to 1.6 m thickness that comprises 4- to 6-m-stacks of three to four cycles (Figure 4). Shorter cycles show one of three motifs: gradual coarsening upwards with sharp tops; sharp basal increases, followed by gradual fining; or symmetrical gradual coarsening then fining. These are all interpreted as parasequences (Van Wagoner et al., 1988; Bohacs et al., 2022a) and stack to form parasequence sets that are overall coarsening (prograding), gradually backstepping or fining (retrograding), or, rather imperceptibly, with no overall change (aggrading) — see Van Wagoner et al. (1990), Bohacs et al. (2022b).

↗

Figure 4. Elemental geochemistry showing calcium, rubidium and zirconium data. Grain-size proxies silicon/aluminium and zirconium/rubidium are also set against the Robin Hood’s Bay graphical log, simplified from Figure 2. Stratal parasequence stacking patterns are labelled with arrows representing prograding or aggrading parasequence sets (green) and retrograding sets (blue). The two proxies for grain size show good agreement with textures recorded in the log that are based on visual description s in the field (for detailed field observations see Hesselbo & Jenkyns,1995; van Buchem & Knox,1998; McIlroy, 2007; Hesselbo et al. ,2020a). Stacking patterns, with missing stacks of retrogradational parasequence sets , indicate the principal hiatuses, and are marked by bold red dashed lines at 15 m and 0 m height. Pale dashed red line s mark the bases of retrogradational parasequence sets , where underlying progradational sets are still present. D ashed blue lines mark surfaces of maximum flooding and correspond to an upward transition from retrogradational parasequence sets to progradational parasequence sets. Hgt=height. For full ammonite chronozone names see Figure 2.

4.2.2. Redox sensitivity

Redox-sensitive elements fluctuate primarily in response to changing redox conditions in the water column or at the sediment-water interface. Molybdenum (Mo) is a commonly used redox-sensitive element that is insoluble in reducing, euxinic conditions (Algeo, 2004; Algeo & Maynard, 2004). Other elements such as nickel (Ni) and chromium (Cr) are indicative of less severe dysoxic conditions (Tribovillard et al., 2006). Molybdenum is below detection limit for the majority of the sectionbut is identified consistently from +10 m (Figure 3). In the lower Pliensbachian Jamesoni Zone, Mo is recorded in and out of detection limit, with values as high as 26 µg/g, but an average of ~5.8 µg/g. Other redox-sensitive elements such as vanadium (V) and nickel also reach their highest values (480 µg/g and 133 µg/g, respectively) in the Jamesoni Zone (compatible with previously published lab XRF data for vanadium; van Buchem et al., 1992).

There is relatively little change in iron abundance (see Supplementary Data File 2) with values ranging 4 wt% to 6 wt% (sampling avoided the common siderite nodules); however, sulphur shows an increase from ~0.1% in Sinemurian values to 1.6% in the Pliensbachian, resulting also in S/Fe ratios that rise to well above 0.15 in the Pyritous Shale. These observations are concordant with previous work on the Lower Jurassic (Toarcian) of the Cleveland Basin (McArthur et al., 2008; Thibault et al., 2018) that support the use of simple total S/ total Fe ratios as a measure of Degree of Pyritization (DOP-T), an index that is sensitive to bottom water oxygenation (Raiswell et al., 1988). Thus, DOP-T indicates much higher pyritization in the Jamesoni Zone and Pyritous Shale Member, compatible with the sedimentological and palaeoecological observations for more oxygen-deficient bottom waters.

4.3. Visual and spectral analysis of the stratigraphic depth series

The results of time series analysis of the 44 m long succession from the mid Sinemurian to early Pliensbachian are presented for those palaeoenvironmental proxies described above that showed promising visual cyclicity (Figures 3–5).

↗

Figure 5. 2πMTM power spectra for nickel, δ¹³C_org, Si/Al, and TOC stratigraphic depth series, post-processing (interpolated at 10 cm [organic matter] and 5 cm [HH-XRF] and detrended [1540 cm LOESS ]) . The power spectra show peaks in relative power at 6–4 m, 1.5 m, 58–50 cm and 30–26 cm across all four datasets illustrated.

The nickel, δ¹³C_org, Si/Al and TOC data, which provide the most continuous time series that also reflect depositional conditions, were first analysed in the stratigraphic depth domain (Figure 5). The results show elevated power at 6–4 m (>99.9% CL), 1.6–1.5 and 1.16 m, and 58–50 and 31–26 cm (>95% CL). The visual analysis, shown in Figure 3, suggests 8 to 8.5 longer period cycles, based on successive minima and bundling in δ¹³C_org and TOC, with five cycles in the 5–6 m range and four in the 4–5 m range; one cycle is incomplete at the base of section. In the spectral analysis, for the ~1.5 m band, the highest power occurs at ~1.46 m (Figure 5).

Bandpass filters were used to extract wavelets of interest from the depth series data for analysis, quality control, and tuning. An elliptic filter was created to encompass 1.6–1.0 m, and a second wide, Gaussian filter centred on ~5 m (6.25–4.16 m) was also extracted, based on the power spectral analysis. These filters show excellent agreement with the data and a strong 1:4 ratio that is consistent with the hypothesis of a long orbital eccentricity/short orbital eccentricity relationship (i.e., 405 kyr and 100 kyr cycles). The filter outputs are shown with the raw data in Figure 6.

↗

Figure 6. Robin Hood’s Bay comparisons between depth-series datasets. Bandpass filters in green represent the hypothesised 405 kyr cycles and red represent 100 kyr cycles, taken from MTM power spectral analysis of respective data series. Visual bundling of cycles matches well the bandpass filters. Red shaded intervals are zones suspected to have missing parasequences. In both cases, the 405 kyr cycles are shorter than average and have fewer than the expected four constituent 100 kyr cycles. This suggests ~100 kyr missing at the Double Band bed, and 100 kyr at the East Scar level. The durations of chronozones are estimated by counting the 405 kyr and 100 kyr cycles, and the range accounts for the amount of potentially missing time, assuming that there are no whole 405 kyr cycles missing.

However, there are places where visual analysis does not match the interpreted 405 kyr cycles from bandpass filters. For example, between -12 m and -15 m and between 1 m and 4 m in the section, there are shorter hypothesised 405 kyr cycles (4–5 m thick) and only three constituent cycles making up the bundles. This suggests that candidate 100 kyr cycles are missing (Figure 6). Coarser-grained sediment and concentrations of thick-shelled bivalve fossils (Figure 2) support the hypothesis of localised sediment bypass.

Time-depth conversion was attempted by cycle-counting using 405 kyr and the 100 kyr cycles with the knowledge that differences may occur due to <405-kyr-scale discontinuities or severe condensation. The proxies Ni and δ¹³C_org were chosen for tuning to provide evidence of forcing in both elemental and organic geochemical proxies (Figure 7).

The 405-kyr tuned Ni power spectrum (Figure 7B) appears to better represent the full range of Milankovitch periodicities, with peaks above the 99.9% CL at 405, 125–95, 41 and 35 kyr. Peaks above the 99 % CL occur at 20.6 kyr to 19.6 kyr and 16 kyr. In addition to these, low period cycles at 1.23 Myr and 2.07 Myr were identified in the tuned and un-detrended datasets. Given the spread of power around 40–43 kyr, uncorrected sedimentation rates are still affecting the age model.

The 100-kyr tuned δ¹³C_org power spectrum (Figure 7C) has sharp peaks at 21 kyr, 104–96 kyr, and a broad peak from 396–530 kyr. The 100 kyr tuning (based upon a 100 kyr filter) indicates that the total section duration is 3250 kyr. This is distributed with 500 kyr in the partly sampled Obtusum, 900 kyr in the Oxynotum, 1350 kyr in the Raricostatum zones, and 500 kyr in the partly sampled Jamesoni zone (Table 1).

The whole series tuned to the 405 kyr cycles has a duration of 3550 kyr. This suggests, compared with the shorter 100-kyr tuned time series model, that 300 kyr is missing from the sedimentary record.

↗

Table 1. Estimated durations for the studied section and ammonite chronozones based on the two age models (tuned using 405 kyr and 100 kyr cycles). Differences occur in the Raricostatum, Oxynotum and Obtusum zones, interpreted here as a result of hiatuses.

Ammonite Chronozone	Short Eccentricity Age Model (kyr)	Long Eccentricity Age Model (kyr)	Difference (kyr)
Jamesoni (partial)	500	500	0
Raricostatum	1350	1450	100
Oxynotum	900	1000	100
Obtusum (partial)	500	600	100
Total section	3250	3550	300

↗

Figure 7. (A) Tuned data time series using the 405 kyr cycles (5.12 m bandpass filter). Filters centred on 405 kyr and 100 kyr from the Ni time-series are displayed alongside the Ni proxy. Dashed grey lines represent ammonite chronozone boundaries; dashed red lines are inferred significant hiatus levels. (B ) 2πMTM power spectrum for Ni data tuned using the 405 kyr cycle. (C) 2πMTM power spectrum for δ¹³C_org data tuned using the 100 kyr cycle.

5. Discussion

5.1. Sequence stratigraphy and relative sea-level change

The trends observed in the grain size proxies above are parasequences that each represent gradual changes in sediment composition through time, commonly bounded by surfaces indicative of minor stratigraphic condensation (Van Wagoner et al., 1988; Bohacs et al., 2022a). The stacking of these parasequences to form larger coarsening or fining upward successions shows progressive and graduated changes in the overall energy of the depositional environment (Van Wagoner et al., 1990; Bohacs et al. 2022b). The stacking pattern also highlights stratigraphic levels where the relatively conformable sequence is broken. A typical depositional sequence cycle in the data shows progressive backstepping or retrogradational parasequence sets, followed by an equally progressive set of aggrading then prograding parasequences (Van Wagoner et al., 1990). Alternatively, and in parallel, the succession can be interpreted in the context of five Transgressive-Regressive (T-R) cycles (see Embry, 1993, for general discussion).

In the lowest complete depositional sequence (Figure 8) from the coarser Lower Triplet bed in the Denotatus Subzone of the Obtusum Zone, the fining upwards, backstepping or retrogradational trend is indicative of a decrease in the overall energy regime, and is followed by aggradation and progradation, culminating in the Double Band bed of the Oxynotum Subzone (see Hesselbo et al., 2020a). The second depositional sequence starts with an abrupt drop to low Si/Al values, with inferred missing retrogradational parasequences, and then a continued progradational parasequence set upwards to the Upper Triplet bed. The third depositional sequence shows an expanded retrogradational parasequence set, indicative of overall decreasing depositional energy followed by equally expanded aggradation and progradation culminating at the East Scar bed (Macdonnelli Subzone). Above East Scar, in the fourth depositional sequence, another stack of retrograding parasequences is inferred to be missing based on the sharp jump back to progradation (Macdonnelli Subzone). Above this, a complete and expanded retrogradational parasequence set continues into the Jamesoni Zone and comprises the start of a fifth depositional sequence. The precise distinction of aggradational parasequence sets from either retrogradational or progradational sets can be somewhat subjective, and in Figure 4 we have chosen not to label aggradational parasequence sets.

The absence or severe reduction of the expected retrogradational parasequence sets indicates that the relatively continuous sequence of deposition is broken, especially at the top of the Double Band (Oxynotum Subzone) and the top of East Scar (Macdonnelli Subzone). These stratigraphic levels are candidates for missing strata with durations below the level detectable by ammonite subzone biostratigraphy. Although relatively short, these intervals potentially affect the interpreted duration at the stage level. The absence of retrogradational parasequence sets indicates that the cause of the hiatus is likely to be relative sediment starvation combined with sediment bypass in a relatively shallow shelf setting during relative sea-level rise and/or transgression.

In the much-studied Dorset (southern UK) succession, there are two major hiatal surfaces in the time-equivalent strata: the beds locally known as the Coinstone and the Hummocky (Hesselbo & Jenkyns, 1995). The Coinstone bed evidences removal of the whole Oxynotum Zone, and the Hummocky bed lies at the Pliensbachian–Sinemurian boundary, with two subzones of the Raricostatum Zone (Macdonnelli and Aplanatum) missing (Lang & Spath, 1926; Lang, 1945; Hallam, 1969; Sellwood, 1972; Hesselbo & Palmer, 1992). Both these beds are characterised by hiatus concretions, with the Coinstone showing clear boring and reworking of septarian concretions, evidence of exposure on the seafloor. The origin of these hiatuses has been proposed to be due to physical erosion at sea-level lowstand (Hallam, 1969), or sediment starvation with bio-erosion at sea-level highstand (Hesselbo & Palmer, 1992; Hesselbo & Jenkyns, 1998).

It has been argued from the visual descriptions of outcrops that the large-scale cycles in relative sea level are broadly synchronous between sedimentary basins in the UK (Hesselbo & Jenkyns, 1998; Hesselbo, 2008). This interpretation agrees with other UK and northern European basins, which all show major flooding events at the lower Oxynotum Zone and in the Jamesoni Zone (Figure 8; Hesselbo & Jenkyns, 1998; Barth et al., 2018; Haq, 2018).

A descriptive approach to sequence stratigraphy is commonly limited by the often-cryptic argillaceous units, which are dominant for large parts of the successions (see Hesselbo & Jenkyns, 1998, with respect to the Lias Group). In using high-resolution elemental geochemistry, insights have been gained into the relative changes within and between fine-grained successions. The analysis here, where grain-size change patterns lie at the border of what is visually obvious at outcrop and what is not, suggests that missing strata in Yorkshire occur repeatedly associated with flooding surfaces and within a framework of Milankovitch scale cyclicity.

↗

Figure 8. Elemental grain-size proxy silicon/aluminium set against the summary graphic log for Robin Hood’s Bay. The stacking patterns shown in Figure 4, and summarised here, are indicative of transgressive-regressive cycles and therefore, trends in relative sea-level, illustrated through blue and green triangles. These trends are previously summarised in Hesselbo & Jenkyns (1998) and Hesselbo (2008). Hiatuses are highlighted by strong red dashed lines at -15 m and 0 m height. Evolutive Harmonic Analysis (EHA) conducted on the Si/Al proxy shows offsets in the major power bands associated with the highlighted surfaces in the sequence stratigraphy, giving confidence in the interpretation of hiatuses ; other sequence stratigraphic surfaces also coincide with abrupt changes in sedimentation rates, and some may also be associated with more minor hiatuses.

5.2. Integration of cyclostratigraphy and sequence stratigraphy

In Hesselbo et al. (2020a), it was hypothesised that parasequences at Robin Hood’s Bay were the expression of the ~100 kyr duration short orbital eccentricity cycle. This hypothesis follows interpretations made by van Buchem & McCave (1989) and van Buchem et al. (1992, 1994), who suggested that subtle lithological variations in the Pliensbachian of the same basin were controlled by Milankovitch cyclicity. The time series and visual analysis of the high-resolution dataset of the present study shows that the dominant cyclicity in the Sinemurian is ~ 1.46 m, shorter than the 2 m from the lower Ironstone Shale Member (van Buchem et al., 1992), suggesting a lower overall sedimentation rate for the Siliceous Shale Member (between 1–1.6 cm / kyr).

As argued above, tuning based upon the hypothesis of short eccentricity (~1.46 m) shows statistically significant peaks in power concentrated in the expected 405 kyr and 21 kyr bands, which is consistent with the interpretation of Milankovitch forcing and confirms the hypothesis of the expression of parasequences at the scale of short eccentricity. The strong eccentricity signal likely originates from rectification of eccentricity modulation of precession by geological processes (Weedon, 2003). The relatively low sedimentation rates in the Robin Hood’s Bay section – hence a thin precession band – and intense bioturbation in all sedimentary facies (Sellwood, 1970; Hesselbo & Jenkyns, 1998; McIlroy, 2007; Powell, 2010), are likely reasons for the reduced power of the precession and obliquity signal in the dataset (see Martinez, 2018). The broad peak that results at the long eccentricity frequency band when tuning to the 100 kyr cycle is indicative of spectral smearing due to changes in the thickness of the cycle through the section (Moore & Thomson, 1991; Weedon, 2003). This reflects the observed missing stratigraphy and varying sedimentation rates documented above. The clearest example of missing parasequences is at the Double Band, which is proposed here to accommodate the majority of this missing time.

Additionally, we note that each of the recognised major flooding surfaces is also associated with shifts in the period of cycles seen in the evolutionary spectrum (Evolutive Harmonic Analysis or EHA) in Figure 8. At the least, this indicates a significant change in the sedimentation rates across these surfaces, and the more severe shifts in cycle frequency might also indicate further minor hiatuses.

The differences in the age models based on 100 kyr and 405 kyr tuning give an estimate of the total missing time across these principal hiatuses of 200–400 kyr. The Double Band interval shows potentially two shortened 405 kyr cycles with only 3 constituent 100 kyr peaks each. This suggests at least 200 kyr of missing time at this level. The missing parasequences above the East Scar level suggest a ~100 kyr hiatus. However, we have low confidence in attributing the amount of missing time to individual hiatuses based on analysis of the Robin Hood’s Bay succession alone, and both the approach of estimating total missing time and where to distribute it is only valid if less than a complete 405 kyr cycle is missing; as a result, the estimates here are a minimum constraint. Further work on high-resolution biostratigraphy and correlation with other sections may help to refine age estimates of missing time (see Weedon et al., 2019; Leu et al., 2024).

5.3. Astronomically forced relative sea level change

The strong imprint of Milankovitch periodicity, particularly short orbital eccentricity, on grain-size proxy datasets in the studied section, is interpreted to be the result of geographically widespread and consistent relative sea-level changes (Hesselbo & Jenkyns, 1998; van Buchem & Knox, 1998), and implies that astronomical cycles were forcing sea level in the Jurassic. Sediment supply can also have a strong impact on sediment stacking patterns and express Milankovitch cyclicity (Van der Zwan, 2002), and climatic effects have been postulated to explain many aspects of Jurassic mudstone lithostratigraphy (Weedon 1986; Huang et al., 2010; Hollaar et al., 2023, 2024). However, the similarities in relative sea level changes across many different basins in the Laurasian Seaway (e.g., Rosales et al., 2006; Hesselbo, 2008; Barth et al., 2018; Francheschi et al., 2022), rule out sediment supply or local subsidence and uplift factors in favour of global sea level change, at least when it comes to so-called 3^rd order (~1 Myr duration) depositional sequences.

Global sea level changes in the Cenozoic oceans have been linked to Milankovitch cyclicity through mechanisms that connect solar insolation to polar continental ice volume (Maslin, 2016; Westerhold et al., 2020). However, the extent to which continental ice can be regarded as a viable mechanism in the Jurassic is debated. Some alternation between relatively cold and relatively warm climate modes has been identified for the Jurassic based upon oxygen isotope and sedimentological evidence, especially for the Pliensbachian, Aalenian–Bathonian, and Tithonian (e.g., Korte & Hesselbo, 2011; Korte et al., 2015; Price et al., 2016), and the Pliensbachian–Toarcian transgression has in particular been suggested to represent a phenomenon of complete deglaciation (e.g., Ruebsam & Schwark, 2021; Nordt et al., 2022; Bodin et al., 2023). Similarly, recent work has demonstrated the existence of continental ice in high austral latitudes in the mid-Cretaceous (Alley et al., 2019), a time like the Jurassic also thought to be relatively warm and for which significant continental ice build-up is debated.

In counterargument, elevated atmospheric CO₂ estimated from soil carbonate carbon-isotopes of 3.5–10x pre-industrial values (Li et al., 2020) and high sea surface temperatures (Price et al., 2016; Robinson et al., 2017), reduce the likelihood that ice was able to accumulate with substantial mass at the poles during the Sinemurian and early Pliensbachian, although this appears less problematic for the Late Pliensbachian where atmospheric CO₂ is thought to be at ~3x pre-industrial values based on the fossil leaf stomatal index proxy (Steinthorsdottir & Vajda, 2015), and oxygen-isotopes of marine macrofossils are relatively heavy (Korte & Hesselbo, 2011).

Nevertheless, large eustatic sea level changes have previously been confidently interpreted for greenhouse periods and associated with Milankovitch cyclicity (Immenhauser, 2005; Huang et al., 2010; Sames et al., 2016; Haq, 2017; Lin et al., 2020). In greenhouse worlds, the effect of ‘limno-eustasy’ (water stored on the continents) has been suggested as the major driver of eustatic sea level variations as an alternative to glacioeustasy (Wagreich et al., 2014; Sames et al., 2016). Correlations between marine and mega-lake deposits for the Jurassic are becoming more refined (e.g., Xu et al., 2017; Franceschi et al., 2022), but they are not yet sufficiently accurate and precise to distinguish synchrony or asynchrony of sea-level and lake-level changes, and thus not able be used to determine the relative contributions of potential mechanisms that force global sea level (cf. Franceschi et al., 2022).

5.4. Astrochronological implications

Using the 405 kyr astronomically tuned age model (Figure 7), the durations of the Oxynotum and Raricostatum ammonite chronozones are 1.00 Myr and 1.45 Myr, respectively. These results may be compared with the much longer but lower resolution astrochronology developed from the Mochras borehole δ¹³C_org curve, Cardigan Bay Basin, Wales (Storm et al., 2020), which show some differences. Due to the constraints on the sample resolution in the Sinemurian at Mochras, making more problematic the identification of obliquity and precession cycles, there is reduced confidence in the identification of the 405 kyr cycles (although higher resolution data has allowed recognition of the shorter cycles in the overlying Pliensbachian; Ruhl et al., 2016; Hollaar et al., 2024). This resolution problem is particularly evident for the Oxynotum Zone at Mochras, where the 405 kyr filter in Mochras matches the data poorly, and where a fault (Woodland 1972; Hesselbo et al., 2013) has likely resulted in missing stratigraphy. Additionally, in the Raricostatum Zone, a change in the thickness of cycles suggests reduced sedimentation rates. The interpretation of the Oxynotum (0.7 Myr) and Raricostatum (1.1 Myr) zonal durations from Mochras are nonetheless of the same general order (Huang, 2018; Storm et al., 2020), albeit shorter than argued for Robin Hood’s Bay.

5.5. Sinemurian–Pliensbachian climate events

The dataset covers two intervals of potential biogeochemical change, the Sinemurian–Pliensbachian Boundary Event and the Liasidium Event. The presented data considerably enhances the stratigraphical characterisation of the GSSP for the base of the Pliensbachian Stage (Meister et al., 2006). The GSSP was selected based on a complete ammonite biostratigraphy with no identifiable hiatuses (Meister et al., 2006). Even though we demonstrate here the location of at least two possible gaps in the Sinemurian, Robin Hood’s Bay is the most biostratigraphically complete section for the Sinemurian–Pliensbachian transition and remains a suitable location to study the carbon-isotope stratigraphy of these two events.

The Sinemurian–Pliensbachian Boundary Event is expressed in Robin Hood’s Bay as a ~4.5 ‰ CIE beginning in the Raricostatum Zone at +3 m height. Correlation with other sections, including the more expanded Mochras borehole (Wales), makes it clear that the excursion is a long-term decrease without an abrupt jump at the Sinemurian–Pliensbachian boundary (Storm et al., 2020). The jump in the δ¹³C_org at +9.5 m in Robin Hood’s Bay coincides with a flooding surface, and the flooding surfaces are shown to be associated with sediment starvation at other levels in the section. This jump is therefore interpreted as a stratigraphic artefact in Robin Hood’s Bay due to condensation, as recently also demonstrated for the Prees core drilled in the Cheshire Basin, ~300 km to the SE (Hesselbo et al., 2023). The presence of the excursion also in δ¹³C_wood and δ¹³C_carb (macrofossil) record (Korte & Hesselbo, 2011; Ullmann et al., 2021) indicates consistent perturbations to the exogenic carbon pool, ruling out local changes in δ¹³C_org as an explanation for the overall excursion but of course, the shape and magnitude of the excursion will also be a response to the mixing of different marine and terrestrial organic components that is evident from pyrolysis data (van Buchem et al., 1992; Suan et al., 2015). In the present study of the Sinemurian–Pliensbachian Boundary Event in the Robin Hood’s Bay section, most of the carbon-isotope excursion coincides with a retrogradational parasequence set (i.e. overall transgression), culminating in a widely recognised maximum flooding surface close to the most negative isotopic values in the Jamesoni Zone (Hesselbo & Jenkyns, 1998; Barth et al., 2018). The Liasidium Event, expressed as a double negative CIE in the Oxynotum Zone at Robin Hood’s Bay (Hesselbo et al., 2020a), likewise coincides with a prominent maximum flooding surface (Hesselbo & Jenkyns, 1998).

The Robin Hood’s Bay section studied here is low in organic matter, even during the onset of the carbon-isotope excursion at the Sinemurian–Pliensbachian boundary. The absence of enrichment in redox sensitive elements compared to average shale, and bioturbated facies through all the upper Sinemurian, indicates oxic conditions. During the Pliensbachian the increase in TOC, degree of pyritization, and redox sensitive elements, allied with sedimentary evidence such as local weak laminations, and subdued benthic communities, points to reducing bottom waters and higher productivity by primary producers. Enrichment in Mo is widely used as an indicator of euxinic conditions (Algeo & Lyons, 2006; Tribovillard et al., 2006). Enrichment with other elements such as V in the absence of Mo is more indicative of sub-dysoxic bottom water (Tribovillard et al., 2006). During the Jamesoni Zone, the absolute values of all these elements are low, with Mo only sporadically present above detection limit. This suggests mostly dysoxic with only sporadic anoxic or euxinic conditions developing during the Jamesoni Zone.

The prolonged isotope excursion at the SPBE has an estimated duration of 1.1 Myr from the onset to the most negative values in Robin Hood’s Bay. This is similar to the duration estimated from other sections (Ruhl et al., 2016; Storm et al., 2020; Pieńkowski et al., 2021; Zhang et al., 2023). High-resolution δ¹³C_org data presented herein show steps that are attributed to the 405 kyr long eccentricity cycle. The recovery from the SPBE is not fully captured in the studied section, but the Mochras borehole indicates a protracted recovery into the Pliensbachian (Ruhl et al., 2016; Storm et al., 2020). There is therefore no evidence of a geologically rapid release of isotopically light carbon, as suggested for other negative CIEs in the Mesozoic (e.g., Hesselbo et al., 2000). The signature instead suggests a steady change in the isotopic composition of the exogenic carbon pools trending to greater proportions of isotopically light ¹²C. Rather than a carbon cycle ‘event’, this appears more like a gradual re-balancing of the carbon cycle that is at least locally climatically modulated. The Liasidium Event CIE is constrained using the astronomical age model to ~ 200 kyr duration and is indicative of more rapid change.

6. Conclusions

Cycles in sediment grain size are expressed as parasequences that are arranged systematically into parasequence sets. The interpretation of parasequence stacking patterns implies hiatuses in the succession, mostly formed during sea-level lowstand and subsequent transgression. The parasequences correspond to short eccentricity (~100 kyr) cycles of astronomical forcing, implying a correlation between relative sea level in the Jurassic and astronomical cycles. The link between eccentricity and relative sea level is interpreted here as being due to eustatic sea-level changes; a conclusion that revives questions about the driving mechanisms for eustasy in greenhouse climates.

Integration of time series analysis with sequence stratigraphy supports the presence of hiatuses, with disruption to regular cyclicity in the intervals of proposed gaps. Analysis suggests the length of the sub-biostratigraphic hiatuses to be less than ~200 kyr in duration, but these estimates are made with low confidence. Most of the missing time at the study site is thought to be in the Macdonelli and Oxynotum subzones. The resultant floating astronomical age model estimates the lengths of the Oxynotum and Raricostatum zones to be 1.00 Myr and 1.45 Myr, respectively, which is ~10% longer than previous estimates.

The late Sinemurian – early Pliensbachian CIE is associated with enhanced organic matter burial and a significant marine transgression. Indicators of palaeo-redox state are compatible with sedimentological and palaeontological observations on the section and suggest that the interval experienced mild dysoxia, with only intermittent anoxic or euxinic conditions, and so better oxygenated than in many other studied sections. Integration with the astronomically tuned age model indicates that the steps observed in the CIE and organic carbon are astronomically influenced, indicating a drawn-out onset of 1.1 Myr from the beginning to the most negative carbon-isotope values. This astronomically paced and relatively long event is not indicative of sudden geological changes to the exogenic carbon pool. Instead, this event appears to be operating on slow timescales indicative of gradual warming and/or changes to the exogenic carbon pool. The Liasidium Event in the Oxynotum Zone also occurs over a period of transgression but does not appear to be associated with dysoxia, at least at the studied location. In contrast to the Sinemurian–Pliensbachian Boundary Event, the Liasidium Event is more rapid, on the order of 200 kyr, estimated from the Robin Hood’s Bay section.

Here we demonstrate that the application of sequence stratigraphy and astrochronology within the same succession, in this case, an argillaceous shallow-marine epicontinental shelf, can provide insight regarding the location of cryptic hiatuses and can result in significant improvements for palaeoclimatic and astrochronological interpretations and age models.

Acknowledgements

AJLH acknowledges PhD studentship funding from the Natural Environment Research Council, NERC GW4+ DTP (NE/L002434/1) and the British Geological Survey (Contract GA/16S/018). SPH and CVU acknowledges NERC grant (NE/N018508/1) Integrated Understanding of the Early Jurassic Earth System and Timescale (JET).

Author contributions

AJLH carried out fieldwork and laboratory analysis and wrote the first draft; CVU provided advice on geochemical analysis and interpretation and edited the manuscript; LAH provided advice on cycle analysis and interpretation and edited the manuscript; KNP revised the ammonite chronostratigraphy; SPH conceived the study, acquired funding, wrote the final draft and provided project supervision.

Data Availability

All relevant data are available in the accompanying supplementary data files.

Conflict of Interest

The authors declare that they have no conflicts of interest in relation the contents of this work.