A new age model for the Pliocene of the southern North Sea basin: a multi-proxy climate reconstruction

The mid-Piacenzian Warm Period (mPWP; 3264–3025 ka) represents the most recent interval in Earth's history where atmospheric CO₂ levels were similar to today. The reconstruction of sea surface temperatures (SSTs) and climate modelling studies has shown that global temperatures were 2–4 ^∘C warmer than present. However, detailed reconstructions of marginal seas and/or coastal zones, linking the coastal and continental climate evolution, are lacking. This is in part due to the absence of precise age models for coastal sedimentary successions, as they are generally formed by dynamic depositional systems with varying sediment and freshwater inputs. Here, we present a multi-proxy record of Pliocene climate change in the coastal southern North Sea basin (SNSB) based on the sedimentary record from borehole Hank, the Netherlands. The marginal marine setting of the Hank borehole during the late Pliocene provides an excellent opportunity to correlate marine and terrestrial signals due to continental sediment input mainly derived from the proto-Rhine–Meuse River. We improve the existing low-resolution palynology-based age model for the Hank borehole using stable oxygen and carbon isotope (δ¹⁸O and δ¹³C) measurements of the endobenthic foraminifera species Cassidulina laevigata, integrated with biochrono- and seismostratigraphy. Identification of hiatuses and freshwater effects in the record allows us to isolate glacial–interglacial climate signals in order to tune the endobenthic oxygen stable isotope record to a global benthic δ¹⁸O stack. This results in a tuned age framework for the SNSB for the late Pliocene (∼3190–2770 ka). Our multi-proxy climate reconstruction for the interval which covers part of the mPWP (∼3190–3000 ka) shows a strong agreement between lipid biomarker and palynology-based terrestrial temperature proxies, which suggest a stable climate, 1–2 ^∘C warmer than present. In the marine realm, however, biomarker-based SSTs show a large range of variation (10 ^∘C). Nevertheless, the fluctuation is comparable to other SST records from the North Atlantic and Nordic Seas, suggesting that a common factor, possibly ocean circulation, exerted a strong influence over SSTs in the North Atlantic and the North Sea at this time.

1 Introduction

The Pliocene Epoch (5.33–2.6 Ma) is a frequently targeted interval for palaeoenvironmental reconstructions because it is considered an analogue for future climate change. For example, atmospheric CO₂ concentrations (380–420 ppmv; Seki et al., 2010; Zhang et al., 2013) and continental configurations during the Pliocene were largely similar to present. Detailed proxy and model comparisons for the so-called mid-Piacenzian Warm Period (mPWP; 3264–3025 ka) have been the focus of the Pliocene Research, Interpretation and Synoptic Mapping (PRISM) group (Dowsett et al., 2010, 2013) and reveal that global temperatures were on average 2–4 ^∘C warmer than present (Haywood et al., 2013). This makes the mPWP an excellent interval to investigate a warmer world associated with the scenarios for our (near) future summarized by the Intergovernmental Panel on Climate Change (IPCC, 2014).

Our understanding of Pliocene climate is largely based on sea surface temperature (SST) reconstructions (e.g. Dowsett et al., 2012), which indicate that global SSTs were 2–6 ^∘C warmer than present. Relatively fewer temperature records exist for the terrestrial realm (Zagwijn, 1992; Salzmann et al., 2013). These records also indicate climate was warmer than present; however, these temperatures are less well constrained due to potential confounding influence of humidity on the temperature reconstructions (Guiot, 1990) and the poorer age control on terrestrial sediment sequences. There are even fewer studies that examine the phase relations and amplitude of variability in coupled land–sea changes (e.g. Kuhlmann et al., 2006), although this information is of key interest for understanding heat transport and the hydrological cycle during the Pliocene. Sediments on continental shelves receive inputs from both the terrestrial and marine environment and would thus contain an archive of land–sea climate evolution. The North Sea basin shelf is a site that potentially hosts a combined record of SST evolution and climate change in the north-western (NW) European continent during the Pliocene due to input of terrestrial material by large European rivers and the active subsidence that provides sediment accumulation space (Gibbard, 1988). Moreover, significant warming of its waters since the second half of the 20th century (0.6 ^∘C rise on average in the period 1962–2001; Perry et al., 2005) indicates the sensitivity of the area for recording climate change. The region has been a type area for Pliocene and early Pleistocene terrestrial stages (see overview in Zagwijn, 1992), but most studied sections lack absolute dating and land–sea correlation, as they target fragmentary fluvial successions (Donders et al., 2007; Kemna and Westerhoff, 2007). However, the shallow marine deposits of the southern North Sea basin (SNSB) allow better chronostratigraphy building through integrated palaeomagnetic, isotope, and biostratigraphic approaches (e.g. Kuhlmann et al., 2006; Noorbergen et al., 2015; Donders et al., 2018).

Despite the promising preconditions that should enable Pliocene climate reconstruction using the sedimentary archive of the SNSB, the generation of an independently calibrated age model for coastal zone sediments is often complicated by complex interactions between sea level, sediment supply, and biotic factors (e.g. Krantz, 1991; Jacobs, 2008; Noorbergen et al., 2015; Donders et al., 2018), which may alter sedimentation rates or cause hiatuses resulting from periods with erosion or non-deposition. Indeed, the Pliocene SNSB was a dynamic system in which multiple westward advances of the Eridanos and Rhine–Meuse rivers generated clinoform successions (Jansen et al., 2004; Kuhlmann and Wong, 2008; Harding, 2015). The sedimentary record thus needs to be critically evaluated on its stratigraphic continuity before it can be compared with records from adjacent areas, such as the Nordic Seas and the North Atlantic. Munsterman (2016) reported a Pliocene-age sequence of coastal marine sediments from Hank, located in the south-west of the Netherlands. The current age framework for the sequence is based on first and last occurrence dates (FODs and LODs, respectively) of dinoflagellate cysts (Dearing Crampton-Flood et al., 2018). Due to the lack of an independent age constraint in the SNSB, FODs and LODs were inferred from those in the Nordic Seas and the North Atlantic, introducing an unknown range of age uncertainty to the biostratigraphic age model (Dearing Crampton-Flood et al., 2018). Furthermore, the resolution of the age model is too low (nine biostratigraphic age tie points for the interval ∼4.5–2.5 Ma) to identify possible hiatuses or changes in deposition, preventing comparison of the record to other archives from the Northern Hemisphere.

The established method for age model construction involves measuring the stable oxygen isotope content (δ¹⁸O) of foraminiferal tests and matching the variability to a global benthic δ¹⁸O reference stack, such as LR04 (Lisiecki and Raymo, 2005). However, in coastal settings this method is complicated due to isotopically lighter freshwater input, which alters the δ¹⁸O value of the foraminiferal tests (Delaygue et al., 2001; Lubinski et al., 2001). Recently, Noorbergen et al. (2015) were successful in creating a tuned age model for the early Quaternary shallow marine interval of borehole Noordwijk, also located in the SNSB, using the δ¹⁸O and δ¹³C values of the endobenthic foraminifera (Bulimina aculeata, Cassidulina laevigata, and Elphidiella hannai). The depth habitat of endobenthic foraminifera in the sediment provides a moderate degree of shelter from disturbances such as reworking by bottom currents and freshwater input. Although vital and microhabitat effects still influenced the absolute δ¹⁸O values of these foraminifera and caused an offset towards more positive values, the trends in δ¹⁸O at Noordwijk clearly resembled those of LR04 (Lisiecki and Raymo, 2005; Noorbergen et al., 2015).

In this study, we follow the approach of Noorbergen et al. (2015) and establish δ¹⁸O and δ¹³C records measured from the endobenthic foraminifera Cassidulina laevigata in the Hank borehole to improve the current low-resolution biostratigraphic age model (Dearing Crampton-Flood et al., 2018) for the Pliocene SNSB. Reconstruction of the age model is further supported by the identification of hiatuses based on seismic information and gamma ray logging. Subsequently, we complement the existing mean annual temperature (MAT) record for Pliocene NW Europe based on soil bacterial membrane lipid distributions stored in the Hank sediments (Dearing Crampton-Flood et al., 2018), with multi-proxy records of SST, prevailing vegetation, and terrestrial input based on lipid biomarker proxies, pollen, and dinoflagellate cysts. This enables us to for the first time directly compare marine and terrestrial climate evolution of the SNSB and continental NW Europe during the mid-Piacenzian Warm Period.

2 Methods

2.1 Geological setting and study site

The Pliocene North Sea was confined by several landmasses except towards the north, where it opened into the Atlantic (Ziegler, 1990). At times, there may have also been a connection via the English Channel to the North Atlantic, indicated by planktonic elements and bryozoan-dominated deposits in the Pliocene Coralline Crag Formation in eastern England (Fig. 1; Funnel, 1996) that bear resemblance to modern deposits. However, the connection of the North Sea to the North Atlantic via the English Channel may only have existed during periods of high sea level (Gibbard and Lewin, 2016). Modern modelled transport estimates from the Hybrid Coordinate Ocean Model (HYCOM) indicate that the total mean inflow at the northern boundary is approximately 14 times higher than that of the English Channel (Winther and Johannessen, 2006). Thus, water inflow/outflow through the English Channel was probably limited, regardless of whether there was an opening to the North Atlantic or not. In addition to a main marine water supply via the North Atlantic, the Eridanos River, draining the Fennoscandian shield, and the proto-Rhine–Meuse River, draining north-western Europe, delivered freshwater to the North Sea in the Pliocene (Fig. 1). The proto-Rhine–Meuse river system existed for a large part of the Pliocene, initially draining the Rhenish Massif and later making a connection with the Alps in the latest Pliocene (Boenigk, 2002). During the Pliocene, the sediment supply by the Eridanos river system to the Ruhr Valley rift system was limited, such that the Rhine–Meuse river system was the predominant supplier of sediments in the study area (Westerhoff, 2009). The water depth of the North Sea during the Pliocene and the Pleistocene was approximately 100–300 m in the central part of the basin (Donders et al., 2018) and approximately 60–100 m in the southern North Sea basin (Overeem et al., 2001).

Figure 1 Pliocene palaeogeography in the North Sea basin (Gibbard and Lewin, 2003; Knox et al., 2010). The location of the Hank borehole is represented by a red star. Major rivers and sediment inputs are represented by blue and orange arrows, respectively. Other locations mentioned in the text are indicated. Figure modified from Gibbard and Lewin (2016).

The study site (51^∘43^′ N, 4^∘55^′ E) is located within the current Rhine–Meuse–Scheldt delta in the region around Hank, the Netherlands. The Hank site is located within the Ruhr Valley Rift, a region that experienced relatively high tectonic subsidence during the late Cenozoic (Van Balen et al., 2000). The current drainage area of the Rhine–Meuse–Scheldt river system is 221 000 km²; however it was likely smaller in the Pliocene (van den Brink et al., 1992; Boenigk, 2002). In 2001, air-lifting well technology was used to drill the Hank borehole (B44E0146) to a base of 404 m. Intervals were drilled every 1 m such that each sample taken from the metre intervals is an integrated mixture. The gamma ray log of the borehole is readily accessible from an online database (https://www.dinoloket.nl, last access: 15 April 2019). In addition, a seismic section is available and covers an east–west transect of the Meuse River (Maas2002 survey, https://www.nlog.nl, last access: 10 February 2019). The lithology of the Hank borehole (Fig. 2a) is described by the Geological Survey of the Netherlands (TNO) and Dearing Crampton-Flood et al. (2018). In short, the base of the succession corresponds to the upper part of the shallow marine Breda Formation, followed by the sandy, occasionally silty, and clay-rich marine delta front deposits sometimes containing shell fragments, so-called “crags”, belonging to the Oosterhout Formation. The overlying Maassluis Formation contains silty shell-bearing deltaic-to-estuarine deposits. For this study, the interval 404–136 m was considered, covering ∼4.5–2.5 Myr based on the biostratigraphic age model of Dearing Crampton-Flood et al. (2018).

Figure 2 Marine proxies for the Hank borehole. (a) The depth and lithology of the Hank sediments, with shell material qualitatively indicated by shell symbols. (b) The smoothed gamma ray (GR) log (https://www.dinoloket.nl, last access: 15 April 2019). (c) Stable oxygen and (d) stable carbon isotope records for the endobenthic Cassidulina laevigata. (e) SST records based on TEX₈₆ (blue diamonds), $U_{\mathrm{37}}^{K^{\prime }}$ (red triangles), and LDI (green squares). (f) Percent cold taxa of dinoflagellate cysts. Intervals 1, 2, and 3 discussed in the text are indicated by green (early Pliocene), grey (mid-Pliocene), and purple (late Pliocene–early Pleistocene). The tuned interval and the position of the hiatus marking the M2 glacial event are represented by a black line.

2.2 Stable isotopes

Deep sea δ¹⁸O and δ¹³C records generally oscillate in antiphase during the Quaternary, due to the waxing and waning of large ice caps on the Northern Hemisphere (e.g. Ruddiman, 2001). During glacial periods, benthic foraminifera incorporate relatively more δ¹⁸O in their calcite test, since more ¹⁶O has been stored in the ice sheets and bottom water temperatures are cold. At the same time, they also incorporate more ¹²C because the total amount of vegetation on land has been reduced during glacial periods, causing an enrichment of the dissolved inorganic carbon (DIC) pool of the oceans. Besides, a reduced thermohaline circulation (THC) during glacial periods may have reduced the contribution of ¹²C-depleted North Atlantic Deep Water to the deep ocean with respect to the ¹²C-enriched Antarctic Bottom Water component. Hence, for the intervals in which the trends in δ¹⁸O and δ¹³C records of the Hank borehole move in antiphase, we interpreted these changes as being the result of glacial–interglacial variability. However, when the δ¹⁸O and δ¹³C are not inversely related, i.e. show a positive correlation, other factors such as riverine freshwater inflow, reworking, and diagenetic influences were most likely the dominant control.

Sediment samples (n=269) from the interval 404–136 m were washed and passed over a series of sieves, after which the > 125 and > 63 µm fractions were collected and dried at 40 ^∘C. Well-preserved (i.e. shiny tests) foraminifera of the endobenthic species Cassidulina laevigata of around the same size were picked from the > 125 µm fraction. Due to the scarcity of foraminifera in some samples, tests were left uncrushed in order to conserve enough material for isotope analysis. The foraminifera were washed ultrasonically in water before weighing. Between 10 and 60 µg of intact tests were weighed per sample. The δ¹⁸O and δ¹³C values were measured on a Thermo Scientific Gas Bench II (Thermo Fisher Scientific) connected to a Delta V mass spectrometer. An in-house NAXOS standard and an internationally accepted NBS-19 standard (δ¹⁸O $=-\mathrm{2.20}$ ‰; δ¹³C =1.95 ‰) were used to calibrate measured isotope ratios to the Vienna Pee Dee Belemnite (VPDB) standard.

2.3 Palynology

Organic-walled dinoflagellates that form a cyst during their life cycle are referred to as dinocysts, and they are preserved in sediments (Head, 1996). Dinocyst assemblages in marine sediments are investigated to infer environmental parameters such as temperature and productivity in surface waters (Rochon et al., 1999; Zonneveld et al., 2013) and can be used as such to reconstruct past climate changes in downcore sediment records (e.g. Pross and Brinkhuis, 2005; Hennissen et al., 2014, 2017). Terrestrial palynomorphs are derived from vegetation and are delivered to coastal marine sediments through wind or river runoff. The pollen and spore (or sporomorph) assemblage in a downcore sediment record like the coastal marine Hank site can thus indicate the type of vegetation in the nearby continent, which can then be used to infer precipitation and/or temperature regimes of the source area (e.g. Heusser and Shackleton, 1979; Donders et al., 2009; Kotthoff et al., 2014).

Standard palynological techniques were used to process 82 selected samples. HCl and HF digestion followed by 15 µm sieving was carried out according to Janssen and Dammers (2008). Both dinocysts and sporomorphs were counted under a light microscope at 400× magnification until a minimum of 200 specimens were found. Rare species were identified during a final scan of the microscope slide. For dinocysts, the taxonomy of Williams et al. (2017) is used.

Some dinocyst taxa prefer cooler (sub)polar waters; hence we take the sum of those taxa and use that as an indicator for SST in the SNSB. We calculate the percent of cold-adapted dinocysts as the sum of the following species, including Bitectatodinium spp., Habibacysta tectata, Filisphaera filifera, Headinium spp., Filisphaera spp., Islandinium spp., Habibacysta spp., Islandinium euaxum, and Bitectatonium tepikiense, over the sum of all dinocysts in the Hank borehole, following the approach adopted by Versteegh and Zonneveld (1994), Donders et al. (2009, 2018), and De Schepper et al. (2011).

A subset of 25 samples was analysed for detailed pollen assemblages to provide independent long-term trends in climate and vegetation cover. Late Neogene pollen types can, in most cases, be related to extant genera and families (e.g. Donders et al., 2009; Larsson et al., 2011). Percent abundances are calculated based on total pollen and spores excluding bisaccate taxa, freshwater algae, and Osmunda spores due to peak abundance in one sample of the latter. Bisaccate pollen abundances are excluded because they are heavily influenced by on- to offshore trends (Mudie and McCarthy, 1994) and therefore do not primarily represent tree abundance.

The terrestrial / marine (T / M) ratio of palynomorphs takes the sum of all sporomorphs and divides by the sum of all sporomorphs and dinocysts, $\mathrm{T}/(\mathrm{T}+\mathrm{M})$. The sum of sporomorphs excludes bisaccate taxa. The T ∕ M ratio is commonly used as a relative measure of sea level variations and therefore distance to the coast (e.g. Donders et al., 2009; Kotthoff et al., 2014).

2.4 Lipid biomarkers and proxies

We use three independent organic temperature proxies for sea surface temperature based on different lipid biomarkers. The TEX₈₆ is a proxy based on the temperature sensitivity of isoprenoidal glycerol dialkyl glycerol tetraethers (isoGDGTs), membrane lipids of marine archaea (Schouten et al., 2002). An increase in the relative abundance of isoGDGTs containing more cyclopentane moieties was found to correlate with SSTs (Schouten et al., 2002). Here we use the global core-top calibration of Kim et al. (2010) to translate TEX₈₆ values to SSTs. Since isoGDGTs are also produced in soils, albeit in minor amounts, they may alter the marine temperature signal during periods with large contributions from land. The relative input of (fluvially discharged) terrestrial organic matter (OM) can be determined using the ratio of branched GDGTs (brGDGTs), which are produced in soils (Sinninghe Damsté et al., 2000; Weijers et al., 2007) and rivers (Zell et al., 2013), with crenarchaeol, an isoGDGT exclusively produced by marine Thaumarchaeota (Sinninghe Damsté et al., 2002). This ratio is quantified in the Branched And Isoprenoid Tetraether (BIT) index (Hopmans et al., 2004), where high BIT indicates a high continental OM input and vice versa. A BIT index > 0.3 is generally used as a cutoff for the validity of TEX₈₆-based SST estimates (Weijers et al., 2006). Secondly, the $U_{\mathrm{37}}^{K^{\prime }}$ index is used as a proxy for SST based on the degree of unsaturation of C₃₇ alkenones produced by marine haptophyte algae (Prahl and Wakeham, 1987). An increased abundance of the tri-unsaturated relative to the di-unsaturated C₃₇ alkenones, expressed as the $U_{\mathrm{37}}^{K^{\prime }}$ index, is linked with decreasing temperature, an adaptation thought to retain membrane fluidity in cooler environments (Marlowe et al., 1984). $U_{\mathrm{37}}^{K^{\prime }}$ values can be converted to SSTs using the global core-top calibration of Müller et al. (1998), with a calibration error of 1.5 ^∘C. Finally, SSTs can be reconstructed based on the relative distribution of long-chain diols, which are dihydroxylated lipids with 22–38 carbon atoms. The C₂₈ 1,13-, C₃₀ 1,13-, and C₃₀ 1,15-diols are most commonly found in seawater and have a putative phytoplankton source (Volkman et al., 1992; Rampen et al., 2007, 2011). The distribution of these three diols is used to formulate the long-chain diol index (LDI), which can be converted to SST using the calibration of Rampen et al. (2012), with a calibration error of 2.0 ^∘C. Furthermore, since freshwater eustigmatophyte algae produce C₃₂ diols (Volkman et al., 1992, 1999; Gelin et al., 1997), the percentage of the C₃₂ diol versus that of the marine C₂₈ 1,13-, C₃₀ 1,13-, and C₃₀ 1,15-diols can be used as an indicator for freshwater discharge (Lattaud et al., 2017).

Lipid biomarkers were previously extracted from the sediments (n=155) and separated into polarity fractions (Dearing Crampton-Flood et al., 2018). The polar fractions, containing GDGTs, were analysed on an Agilent 1260 Infinity ultra-high performance liquid chromatography instrument (UHPLC) coupled to an Agilent 6130 single quadrapole mass detector with settings following Hopmans et al. (2016). For more details on method and solvent programme see Dearing Crampton-Flood et al. (2018). Selected ion monitoring (SIM) was used to detect [M-H]⁺ ions of the isoGDGTS: m∕z 1302, 1300, 1298, 1296, 1292.

After GDGT analysis, polar fractions were silylated by the addition of N,O-bis(trimethylsilyl)trifluoroacetamide (BSTFA) and pyridine (60 ^∘C; 20 min). A Thermo Scientific TRACE gas chromatograph (GC) coupled to a Thermo Scientific DSQ mass spectrometer (MS) was used to analyse long-chain diol distributions in SIM mode (m∕z 299, 313, 327, and 341) at the NIOZ Royal Netherlands Institute. The temperature programme was as follows: 70 ^∘C for 1 min, then ramped to 130 ^∘C at 20 ^∘C min⁻¹, then ramped to 320 ^∘C at 4 ^∘C min⁻¹, and then held for 25 min.

Ketone fractions, containing the C₃₇ alkenones, were analysed using gas chromatography with flame ionization detection (GC-FID). Samples were injected (1 µL) manually on a Hewlett Packard 6890 series GC system equipped with a CP-Sil 5 fused silica capillary column (25 m×0.32 mm; film thickness 0.12 µm) and a 0.53 mm precolumn. The oven temperature programme was similar to that used for long-chain diol analysis.

3 Results

3.1 Stable isotopes of Cassidulina laevigata

Foraminifera preservation in the intervals 404–386 and 204–136 m was either very low or nonexistent. Furthermore, the low abundance and poor preservation of foraminifera in the crag material (220–205 m in Fig. 2c) also presented a challenge for picking. The δ¹⁸O_cass. values of foraminifera in the sediment from the remaining intervals (n= 136) range from −1.0 ‰ to 3.6 ‰ (Fig. 2c). The variability in δ¹⁸O_cass. values between maxima and adjacent minima in the Hank record ranges from ∼1 ‰ to 4 ‰. The stable carbon isotope record varies between −3.8 ‰ to 0.6 ‰ (Fig. 2d), and the variability is ∼0.3 ‰ to 2.3 ‰. Discounting the sample at 206 m, the variability in the δ¹³C_cass. record is approximately ∼1 ‰ (Fig. 2d). The δ¹⁸O_cass. and δ¹³C_cass. records show antiphase relationships in the intervals 378–360, 345–310, 304–292, 282–279, 275–270, 252–242, 232–227, and 224–218 m (Fig. S1a in the Supplement).

3.2 Seismic profile

The ∼15 km east to west seismic profile of the Meuse River, including the location of the Hank borehole, spans a depth of > 500 m (Maas2002 survey, https://www.nlog.nl, last access: 10 February 2019; Fig. 3). Comparison of the formations of the Hank borehole with the seismic depth profile in Fig. 4 indicates that the Breda Formation at 404–370 m is characterized by horizontal-reflection patterns, likely indicating shallow marine conditions. The eastern continuation of the seismic line reveals that these horizontal strata can be interpreted as shelf toesets of westward-prograding deltaic clinoforms.

Figure 3 Seismic east to west depth profile at the Meuse River (from Maas2002 survey, https://www.nlog.nl, last access: 10 February 2019) with the location of the Hank borehole (B44E0146) and corresponding formations indicated. The smoothed gamma ray log (from https://www.dinoloket.nl, last access: 15 April 2019; white) and lithology of the borehole are provided for context. Stars and age ranges refer to the biostratigraphic age model of Dearing Crampton-Flood et al. (2018). The orange and yellow lines represent the boundaries of the Breda and Oosterhout (revised) and the Oosterhout and Maassluis formations, respectively.

Figure 4 Terrestrial proxies for the Hank borehole. (a) The depth and lithology of the Hank sediments, with shell material qualitatively indicated by shell symbols. (b) The smoothed gamma ray (GR) log (from https://www.dinoloket.nl, last access: 15 April 2019) and (c) the relative input of terrestrial organic material to the Hank sediments based on the terrestrial / marine ratio of palynomorphs (black line), the Branched and Isoprenoid Tetraether (BIT) index (orange circles), and the percent of C₃₂ diol (blue squares). Pollen records expressed as percent of total pollen: (d) other conifers (brown line), Taxodioidae (orange line) and MAT (in degrees Celsius; green line; Dearing Crampton-Flood et al., 2018), (e) angiosperm trees, (f) herbs, (g) heather, (h) spores, (i) Pinus, and (j) Osmunda. Intervals 1, 2, and 3 discussed in the text are indicated by green (early Pliocene), grey (mid-Pliocene), and purple (late Pliocene–early Pleistocene). The tuned interval and the position of the hiatus marking the M2 glacial event are represented by a black line.

3.3 Palynology

The palynomorphs in the Hank sediments are well preserved. The borehole can be divided into three main intervals according to the (co)dominance of the marine/terrestrial palynomorphs: 1, 2, and 3, which roughly correspond to the (1) early Pliocene, (2) mid-Pliocene, and (3) late Pliocene–early Pleistocene. In the deepest part of the borehole (404–330 m; interval 1), the marine component of the palynomorph assemblage clearly exceeds the terrestrial, as evidenced by the low T ∕ M values (Fig. 4c). An isolated sporomorph peak and (sub)polar dinocyst peak are visible at 383 m (Figs. 2f, 4c). Interval 2 from 330 to 187 m shows a fluctuating ratio between the marine and terrestrial elements (Fig. 4c). The cold-adapted dinocysts also show fluctuations indicating alternating warmer and cooler periods (Fig. 2f). One striking feature is the increase in cold-adapted dinocyst abundance and simultaneous Osmunda acme at 305 m (Figs. 2f, 4j). The third, interval 3, spans the upper part of the borehole (187–136 m), and sporomorphs in particular dominate this interval, visible by the consistently high T∕M (Fig. 4c). Interval 3 shows an increased occurrence of coastal marine genera, like Lingulodinium. The increased gamma ray values at ∼175 m (Fig. 2b) are the result of the abundance of (shell) concretions, not clays, and as such, do not indicate a more distal environment but rather a development toward a more proximal environment. At 154–153 m, the marine indicators in the borehole are reduced to just 0.5 % of the total sum of palynomorphs. However, the highest abundance of cold-adapted dinocysts, mostly composed of taxa like Habibacysta tectata, is at 154 m (Fig. 2f). This depth is also marked by the complete disappearance of the dinocyst genus Barssidinium spp. with an LOD at 157 m (Dearing Crampton-Flood et al., 2018). (Sub)tropical species like Lingulodinium machaeophorum, Operculodinium israelianum, Spiniferites mirabilis, Tectatodinium pellitum, and Tuberculodinium vancampoae are missing at this depth. The uppermost (154–136 m) interval indicates an estuarine to deltaic environment, due to the presence of freshwater and brackish water algae species Pediastrum and Botryococcus. In contrast, the freshwater indicators are (almost) absent in intervals 1 and 2. The assemblages of interval 3 are also characterized by a fluctuating abundance of cold-adapted dinocysts (Fig. 2f).

The pollen assemblages are dominated by tree pollen, particularly conifers – Pinus, Picea, Abies, Taxodioidae-type (including Glyptostrobus and Taxodium), Sciatopitys, and Tsuga – but with increasing proportions of grasses (Poaceae, Cyperaceae) and heather (Ericales) in interval 2 and significantly increased amounts of fern spores from 260 m and up (Fig. S1). Exclusion of the bisaccate types from the percentage sum (to counter effects of sea level change on the diverse transport capacity of pollen; Neves effect) causes the percentage of pollen sum to be relatively low (∼100 grains/sample). The current pollen record should therefore be used for main quantitative trends and not to delineate individual ranges of taxa. The angiosperm tree abundances average about 20 % and shows no significant long-term change towards the top of the sequence (Fig. 4e). The angiosperm tree pollen record is diverse, although few taxa are continuously present and consists mostly of Quercus robur-type with significant proportions of Pterocarya, Fagus, Carpinus and, above 240 m, Ulmus (Fig. S1). The Taxodioidae-type shows a distinct long-term decline superimposed by three shorter minima in the end of interval 2 and beginning of 3, occurring at 205, 235 and 170 m (Figs. 4d, S2).

3.4 Lipid biomarkers

IsoGDGTs are present in high abundances throughout the borehole, as evidenced by the high total organic carbon (TOC)-normalized concentrations of crenarchaeol (0.2–130 µg g⁻¹ TOC; Dearing Crampton-Flood et al., 2018). IsoGDGT-based SSTs are reconstructed for those sediments where BIT < 0.3, i.e. between 404 and 219 m (n=66). TEX${}_{\mathrm{86}}^H$-reconstructed SSTs range between 7 and 13 ^∘C but do not show a clear trend over time (Fig. 2e).

Alkenones are present in the majority (n=111) of the samples. However, they are below the detection limit in many of the sediments from 250 to 200 m. Alkenones re-emerge in the interval from 197 to 178 m. The $U_{\mathrm{37}}^{K^{\prime }}$ index values range between 0.30 and 0.83 and correspond to a SST range of 8–24 ^∘C (Fig. 2e). In the early Pliocene (1), the $U_{\mathrm{37}}^{K^{\prime }}$ SST record shows the largest fluctuations in temperature (ΔT=15 ^∘C) and an average SST of 19 ^∘C. Similar variability (ΔT=14 ^∘C) is observed in interval 2 (middle–late Pliocene), although the average SST drops slightly to 17 ^∘C. Alkenones are present around or below the detection limit in interval 3 (late Pliocene–early Pleistocene), so no SSTs can be calculated. Notably, $U_{\mathrm{37}}^{K^{\prime }}$ SSTs show a warming of 11 ^∘C during the interval between 290 and 260 m (Fig. 2e).

Long-chain diols are below the detection limit in a large proportion of the Hank borehole. SSTs can be reconstructed for a select few samples in intervals 1 and 3, and they show scattered temperatures in a range of 13 ^∘C (Fig. 2e). The sediments in interval 2 contain enough diols to enable a semi-continuous SST reconstruction. The range of LDI SSTs in this section is 4–18 ^∘C (Fig. 2e). The record shows a strong warming trend of ∼12 ^∘C from 295 to 263 m, coeval with the trend in the $U_{\mathrm{37}}^{K^{\prime }}$ record (Fig. 2). The percent of C₃₂ diol is generally high (∼36 %–57 %) in intervals 1 and 2 (early–middle Pliocene; Fig. 4c), indicating a modest to strong freshwater input (cf. Lattaud et al., 2017). The percent of C₃₂ diol slightly decreases (10%) over 294–264 m, indicating a gradually decreasing influence of riverine OM and/or an increase in the abundance of the marine C₂₈ 1,13-, C₃₀ 1,13-, and C₃₀ 1,15-diols. In interval 3, the percent of C₃₂ diol exhibits a strong increasing (21–59%) trend (187–136 m; Fig. 4c).

4 Discussion

4.1 Depositional record and unconformities

The changes in depositional environment in the Hank borehole from open marine to coastal marine and successively estuarine conditions are based upon the BIT, TOC, and δ¹³C_org records presented in Dearing Crampton-Flood et al. (2018) and the biological changes in the abundances of typical marine (dinocysts and test linings of foraminifera), estuarine/freshwater algae species, and sporomorph assemblages (Munsterman, 2016), summarized in the T ∕ M ratio. A transition of marine OM during the Pliocene to more terrestrial OM input towards the Pleistocene starts approximately at 190 m as evidenced by these indicators (Fig. 4c). This reflects the increasing influence of the Rhine–Meuse River to the site. This progradation of the Rhine–Meuse River is confirmed by the westward progradation of the depositional (delta) system that can be seen in the seismic section (Sect. 3.2, Fig. 3).

Toward the long-term shallowing trend in the borehole, several depositional changes and unconformable surfaces are recognized in the seismic profile, which need to be taken into account regarding the construction of a valid age model. The transition from the Breda Formation to the overlying Oosterhout Formation is marked by a distinct angular unconformity referred to as the late Miocene unconformity (LMU; Munsterman et al., 2019; Fig. 3). The seismic data of the overlying Oosterhout Formation (middle–late Pliocene) indicate a twofold subdivision. The lower unit of the subdivision at 370–314 m is characterized by convex-downward reflection patterns that correspond to an open marine signature (with corresponding low sedimentation rates). This is confirmed by the transition to finer-grained sediments (silts) over 381–352 m attributed to a more distal setting (Fig. 2a). This interval is also characterized by an increased abundance of dinocysts with a preference for open marine conditions, like the genus Spiniferites. The facies of interval 352–338 m indicates shallow to open marine conditions with a temperate to (sub)tropical SST. A coarsening upward trend is corroborated by the gradual decrease in gamma ray values (Fig. 2b). In the second Oosterhout unit from 314 to 157 m the environment is shallow marine, and several stacked clinoform sets are visible in the seismic profile (Fig. 3). The transition between the two Oosterhout units is clearly visible as a downlap surface around 314 m that corresponds to a sequence boundary and possible hiatus (Fig. 3). This is coupled with a dramatic decrease in sedimentation rate in the initial age model, which places this interval within the scope of the M2 glacial event (∼3.3 Ma; Dearing Crampton-Flood, 2018).

Above the sequence boundary at 314 m, an increase in water depth to ∼80–100 m can be deduced from the height of the clinoforms where the topsets represent the fluvial system debouching at the coastline. This water depth estimate is comparable to the estimate for Pliocene water depth from an integrated seismo-stratigraphic study of the SNSB (Overeem et al., 2001). The topset beds in the Hank borehole show an abundance of shell crag facies material corroborating the coastal setting (Figs. 2a, 3). The stratigraphic stacking is first purely progradational (clinoforms) and changes to both progradational and aggradational higher up, suggesting progressive fluvial influence replacing the marine environment. In the Hank borehole, this change is marked by a distinct clay layer at 292–271 m (Fig. 2a). In the upper part of the second Oosterhout Formation at depths of 260 m and upwards, the increased proportion of heather and grasses is generally considered indicative of colder and drier terrestrial climate (Faegri et al., 1989; Fig. 4). The decrease in Taxodioidae-type pollen over interval 2 (Fig. 4d) and the Oosterhout to the Maassluis formations further indicates a cooling terrestrial climate and is classically recognized as the top Pliocene (late Reuverian) in the continental zonation of Zagwijn (1960), although in the onshore type area the sequences are most likely fragmented time intervals bounded by several hiatuses (Donders et al., 2007).

At the transition of the Oosterhout Formation to the Maassluis Formation, concave-downward reflection patterns may reflect channel incisions into the topsets of the Oosterhout Formation (Fig. 3). The Maassluis Formation (late Pliocene–early Pleistocene; < 2.6 Ma) is composed of horizontal and channel-like strata in the seismic profile (Fig. 3). The environment of the Maassluis Formation becomes more fluvio-deltaic, characterized by the decreased abundance of dinoflagellate cysts and steep rise in the number of sporomorphs, manifested by the high T∕M values (Fig. 4c). Further warm–temperate trees in the Maassluis interval of the record such as Carya, Liquidambar, and Nyssa disappeared in NW Europe in the earliest Pleistocene (Donders et al., 2007), most likely slightly above the level of the top of the Hank sequence. In contrast, land-based studies in the Dutch–German border area (cf. Donders et al., 2007; Westerhoff, 2009) are characterized by relatively more abrupt last occurrences of warm–temperate taxa due to the probably incomplete preservation of the early Pleistocene sequences.

The Plio-Pleistocene transition (2.6 Ma) occurs between 200 and 154 m (Dearing Crampton-Flood et al., 2018). This transition is accompanied by a peak in gamma ray values at ∼175 m (Fig. 2b). However, the coastal marine depositional setting for the Hank borehole in the interval (upper Oosterhout Formation) during the late Pliocene/early Pleistocene above 200 m strongly indicates that these successions are likely not continuous but consist of successive fragments of sedimentation representing short time windows that are bounded by hiatuses (cf. Donders et al., 2007).

Figure 5 Age framework for the Pliocene southern North Sea basin. Correlations of the interglacials of the (a) LR04 stack (Lisiecki and Raymo, 2005) to the interglacials in the (b) δ¹⁸O_cass. record for Hank. (c) Smoothed gamma ray (GR) and (d) lithology, and depth of Hank sediments. Dotted grey lines indicate tie points based on δ¹⁸O values, whereas the dashed line at ∼314 m is based on correlation of sequence boundary with the Poederlee Formation (De Schepper et al., 2009a; Louwye and De Schepper, 2010). Biostratigraphic age estimates (Dearing Crampton-Flood et al., 2018) are shown as black diamonds.

4.2 M2 event

There is global evidence for a large sea level drawdown during the M2 glacial event, but the estimates of the magnitude of its extent vary greatly (38–65 m; Dwyer and Chandler, 2009; Naish and Wilson, 2009; Miller et al., 2011, 2012). However, large uncertainties in the estimation of ice volume prohibit any meaningful estimates of sea level for the Pliocene using the stable isotope measurements of foraminifera (Raymo et al., 2018).

The hiatus at Hank representing the M2 glacial event is also recognized in sequences from the Coralline Crag in the English North Sea (Williams et al., 2009), the Poederlee and Lillo formations in the Belgian North Sea (De Schepper et al., 2009a; Louwye and De Schepper, 2010) and possibly the Nordic Seas (Fig. 1, Risebrobakken et al., 2016). In the Poederlee and Lillo formations of neighbouring Belgium (located 80 km from Hank in the present day), De Schepper et al. (2009) and Louwye and De Schepper (2010) hypothesize that the MIS M2 is correlated with a sequence boundary PIA1 at approximately 3.21 Ma. In addition, equivocal temperature/assemblage signals in the Coralline Crag Formation are hypothesized to be a result of sea-level change associated with the M2 glacial event, which would have decreased or ceased sedimentation entirely (Williams et al., 2009). A Pliocene benthic δ¹⁸O record adjacent to the North Sea (NS) in the Nordic Seas (ODP hole 642B; Risebrobakken et al., 2016) also does not record any strong evidence of the M2 event, and the authors postulated that the M2 event might have occurred during a hiatus in the borehole or may have been a less extreme event in this region compared to other regions (Lisiecki and Raymo, 2005). Due to the M2 event being a globally recognized event (De Schepper et al., 2014), the records from the East of England, Belgium, and the Hank site indicate that a hiatus likely exists over the most acute part of the glacial event in SNSB sediment successions. Thus, the coolest interval (with the presumed lowest sea level and a high hinterland sediment supply) of the M2 event might have not been recorded at Hank because of erosion.

In the sediments occurring above the hiatus marked by the sequence boundary in Fig. 3, large variability in δ¹⁸O_cass. indicates fluctuating climate conditions that may be associated with the onset or the recovery of the M2 event. The fluctuations match those in the records of the BIT index and δ¹³C of organic matter (see Dearing Crampton-Flood et al., 2018), which indicate a closer proximity of the coast to the site, likely as a result of sea level change. The major peak of Osmunda spores (outside of pollen percentage sum) after the hiatus at ∼3210 ka (306 m) could then represent a pioneer phase of marsh vegetation related to a rapid sea level lowering. The (sub)polar dinocyst acme and increase in Operculodinium centrocarpum (Figs. 2, 6a) at 305 m may then represent the restoration of the location of the Hank site to a more distal marine setting within the confinement of the Upper Rhine Graben. The sea level drop associated with the M2 event may have decreased the inflow of Atlantic bottom water currents originating from the northern opening of the North Sea (Kuhlmann et al., 2006). After the M2 event, isostasy may have then strengthened the connection to the North Atlantic (possibly also via the English Channel), which would have allowed the inflow of relatively warmer and saline Atlantic Water fed by the North Atlantic Current (NAC) into the North Sea. The dinocysts of Operculodinium centrocarpum are generally used as a tracer for the NAC in the North Atlantic (De Schepper et al., 2009b; Fig. 6a) and may tentatively be linked to the increasing influence of the North Atlantic to the study site after the hiatus associated with the M2 event. However, Operculodinium centrocarpum is a cosmopolitan species, and in the modern day with a connection via the English Channel to the North Sea, it is not commonly found (Marret and Zonneveld, 2003; Zonneveld et al., 2013). Thus, the re-emergence of Operculodinium centrocarpum over the mPWP interval (Fig. 6) more likely reflects the restoration of marine conditions in the shallow SNSB after the M2 event.

The acme of Osmunda spores coincides with the occurrence of dinocysts characteristic of (sub)polar water masses at ∼3210 ka, further indicating cold conditions (Fig. 6). In addition, the distinct decrease in Taxodioidae-type pollen at the same time indicates that climate conditions were also cold(er) on the continent (Fig. 6c), which is supported by low terrestrial mean air temperatures of 6 ^∘C, independently reconstructed based on brGDGTs (Fig. 6f; Dearing Crampton-Flood et al., 2018). In contrast, all SST reconstructions remain stable during this M2 deglaciation/recovery period (Fig. 6d), suggesting that cold periods on land are better recorded in the sedimentary record than those in the marine realm. Indeed, terrestrial proxies should represent an integrated signal over longer time and larger space (NW Europe), compared to that of the marine proxies, which are confined to the shallow SNSB basin and potentially mostly record warm periods (Sect. 4.3).

Figure 6 Climate proxy records for the southern North Sea basin for the late Pliocene. Age tie points (stars) based on oxygen isotope stratigraphy (black), sequence boundary correlation (red), and biostratigraphy (blue) are indicated. The depth interval covered by the age tying points is 206–330 m. (a) The relative abundance of Operculodinium centrocarpum expressed as a percent total of dinocysts; (b) percent of cold dinocysts; (c) pollen abundances for Taxodioidae (orange) and Osmunda (dark yellow) as a percentage of the pollen sum; (d) SST records based on the TEX₈₆, $U_{\mathrm{37}}^{K^{\prime }}$, and LDI proxies, together with estimates (dotted lines) from oxygen isotope measurements of bivalves (Vignols et al., 2019); (e) the relative input of terrestrial organic material to the Hank sediments based on the Branched and Isoprenoid Tetraether (BIT) index (Dearing Crampton-Flood et al., 2018) and the input of fresh water based on the percent of C₃₂ diol; (f) mean air temperature based on brGDGT-palaeothermometry (Dearing Crampton-Flood et al., 2018); and (g) the benthic oxygen isotope stack of Lisiecki and Raymo (2005).

4.3 Glacial–interglacial variability and tuning

Based upon the age model of Dearing Crampton-Flood et al. (2018), it is clear that the sample resolution is too low to resolve a stable isotope tuning on Milankovitch timescales for the older succession including the Breda and lower Oosterhout formations (404–330 m), but it is sufficient (i.e. < 6 kyr) to resolve individual cycles, in particular above 314 m. Notwithstanding the low resolution between 390 and 314 m, the δ¹⁸O_cass. and δ¹³C_cass. records move in opposite directions, suggesting that global glacial–interglacial ice volume changes were largely influencing the open (shallow) marine environment at Hank during that time. In the upper part of the Hank borehole (314–200 m), the trends in δ¹⁸O_cass. and δ¹³C_cass. do not show a continuous inversed relationship (Fig. S1a), indicating that the ice volume influence is likely obscured by other factors such as riverine freshwater inflow.

The absolute values of the oxygen isotope measurements on Cassidulina laevigata recorded in Hank are substantially lower by approximately 1 ‰–1.5 ‰ than the composite benthic δ¹⁸O values in the LR04 stack (Lisiecki and Raymo, 2005), as well as those of a nearby Pliocene benthic oxygen isotope record from the Nordic Seas (∼2 ‰–3 ‰; Risebrobakken et al., 2016). The offset in absolute values is unlikely due to a species-dependant effect, as δ¹⁸O_cass. values in a nearby Quaternary-age core from Noordwijk (Noorbergen et al., 2015) were comparable to the LR04 stack (Lisiecki and Raymo, 2005). Hence, the relatively low δ¹⁸O_cass. values in the Hank sediments likely reflect the influence of freshwater input at this site, which is proximal to the mouth of the palaeo-Rhine–Meuse River (e.g. Delaygue et al., 2001; Lubinski et al., 2001; Westerhoff, 2009; Fig. 1). Similarly, low benthic and planktic δ¹⁸O values in sediments from the Ionian Sea coinciding with sapropel deposition were attributed to increased freshwater influence at this time (Schmiedl et al., 1998). Furthermore, the large δ¹⁸O_cass. variability in the Hank record (∼1 ‰–4 ‰) compared to that in the LR04 stack record (0.2 ‰–0.7 ‰ during the Pliocene; Lisiecki and Raymo, 2005) indicates that the δ¹⁸O isotopic signature of shallow seawater at the Hank site is sensitive to freshwater input. At the delta front, wave action and winnowing contribute to the mixing of freshwater input in the relatively shallow water column, which explains why an endobenthic species may be affected by freshwater influence, which results in relatively lower absolute δ¹⁸O_cass. values at the Hank site. Thus, salinity changes and sensitivity to freshwater input affect the oxygen isotopes incorporated into Cassidulina species at the Hank site, regardless of the endobenthic habitat.

The upper (Plio-Pleistocene transition) and lower (M2 glacial) stratigraphic boundaries identified in Sects. 4.1 and 4.2 provide a contextual framework to construct a higher-resolution age model for the mPWP (3264–3025 ka) using stable isotopes of Cassidulina laevigata. The open marine signature and relatively horizontally deposited clinoform sets in the second unit of the Oosterhout subdivision from ∼305 to 200 m (Fig. 3) represent a relatively continuous sedimentary record that may be suitable for age model reconstruction, keeping in mind the potential freshwater influence on the δ¹⁸O_cass. record.

In order to use the δ¹⁸O_cass. record for tuning purposes, an understanding of the North Sea hydrogeography and circulation patterns during the Pliocene must be taken into consideration. During cold periods, the North Sea circulation slows due to the reduced sea level and inflow of Atlantic water (Kuhlmann et al., 2006). Stratification in the North Sea due to freshwater input from rivers combined with the sluggish circulation and weak influence of the Atlantic waters make cooler periods problematic to tune to, due to a δ¹⁸O signature that is probably highly localized and erratic. Moreover, Donders et al. (2007) noted that the coldest phase of glacials of the Plio-Pleistocene climate development of coastal areas in the NS is likely to be marked by substantial hiatuses caused by non-deposition and erosion, which may also preclude the use of the transition between the warmer and cooler periods as a tuning anchor. During warmer periods, an increased freshwater input from river outflows is also expected, due to the supposedly wetter climate conditions during interglacials. However, Kuhlmann et al. (2006) linked warmer periods in the Pliocene in the central section of the southern North Sea with the occurrence of Cassidulina laevigata, whose habitat in the modern North Sea is located in the northern part with a strong connection to the Atlantic (Murray, 1991). Thus, tuning the warmer periods in the δ¹⁸O_cass. record at the Hank site with warm periods in the LR04 benthic stack is preferable due to the strong(er) connection to the Atlantic (Kuhlman, 2004), resulting in a relatively more regional signature of the δ¹⁸O_cass. values (Kuhlmann et al., 2006). Moreover, the chance of disturbance/hiatuses that affect the continuity of the sediment record at Hank is decreased in warmer periods.

Using the above reasoning, the sample with the lowest δ¹⁸O_cass. value in each cycle between ∼314 and 200 m in the Hank record can be tuned to the lowest δ¹⁸O value between the M2 event associated with the hiatus at 314 m (∼3.3 Ma) and the Plio-Pleistocene boundary (Sects. 4.1, 4.2) in the LR04 stack, assuming that the low δ¹⁸O values in the Hank borehole represent the warmest part of each interglacial period (Fig. 5). Further investigation into the variation in the δ¹⁸O_cass. cycles in the Hank borehole isotope record reveals unique sawtooth structures, differing from the more symmetrical pattern of cyclicity that is seen in the Pleistocene interval of the LR04 stack. Specifically, cycles G19, G17, and G15 display these reversed sawtooth patterns in the global benthic stack and help pinpoint corresponding cycles in the Hank borehole record (Fig. 5). The reconstructed time window spans ∼3190–2770 ka and thus most of the mPWP. Based on the tuned oxygen isotope age model, the LOD of Invertocysta lacrymosa and Operculodinium? eirikianum in the SNSB can be constrained to ∼3040 and < 2768 ka, respectively (see Dearing Crampton-Flood et al., 2018).

4.4 Late Pliocene climate reconstruction

4.4.1 Marine proxy interpretation

Despite the fact that all three lipid biomarker proxies (TEX₈₆, $U_{\mathrm{37}}^{K^{\prime }}$, and LDI) are calibrated to SST, the records that they generate show remarkable differences and are offset in temperature (Fig. 2e). The Pliocene TEX${}_{\mathrm{86}}^H$ SSTs are 10 ^∘C on average, which is the same temperature as the modern mean SST of the North Sea (Locarnini et al., 2013) and contrasts with other North Sea Pliocene temperature estimates based on ostracod, mollusc, foraminiferal, and dinocyst assemblages (Wood et al., 1993; Kuhlman et al., 2006; Johnson et al., 2009; Williams et al., 2009), all suggesting that the SST of the North Sea was 2–4 ^∘C warmer than present at that time. However, present-day TEX₈₆ reconstructions for core-top sediments in the North Sea range between 4.1 and 9.1 ^∘C (Kim et al., 2010) and thus underestimate the observed modern SST. Lower-than-expected TEX₈₆ values found elsewhere have been explained by a contribution of isoGDGTs produced by a subsurface community (Huguet et al., 2007), but given the shallow water depth (80–100 m) of the SNSB in the Pliocene (Hodgson and Funnel, 1987; Long and Zalasiewicz, 2011; Overeem et al., 2001; this study), it seems unlikely that such a community would have played a role here. This can be further confirmed by calculating the ratio of isoGDGT-2 / isoGDGT-3 ([2] ∕ [3]; Taylor et al., 2013), whose value increases with increasing isoGDGT input from subsurface-dwelling archaea. The [2] ∕ [3] ratio in the Hank borehole is 2.1 on average and always well below the value associated with a deep-water archaea community overprint (> 5; Taylor et al., 2013). Instead, the low TEX${}_{\mathrm{86}}^H$ SSTs are likely a result of seasonal production of isoGDGTs. In the modern North Sea, the main period of thaumarchaeotal blooms and associated isoGDGT production are in the winter months where ammonia is available and competition with phytoplankton is minimal (Herfort et al., 2006; Pitcher et al., 2011), which likely introduces a cold bias in TEX₈₆-based SST estimates for the SNSB. If the TEX₈₆-derived SSTs are interpreted as a winter signal as we argue, then the Hank Pliocene SSTs are approximately 3–6 ^∘C warmer than modern winter SSTs (van Aken, 2008).

Conversely, $U_{\mathrm{37}}^{K^{\prime }}$-reconstructed SSTs are 16 ^∘C on average and thus 2–4 ^∘C higher than the temperature estimates based on ostracod, mollusc, and foraminiferal assemblages (Wood et al., 1993; Kuhlmann et al., 2006; Johnson et al., 2009; Williams et al., 2009) and ca. 6 ^∘C higher than modern annual mean SST. These higher-than-expected $U_{\mathrm{37}}^{K^{\prime }}$ SSTs could in part be caused by a species effect as a result of a contribution from alkenones produced by freshwater haptophyte algae that have little to no correlation of $U_{\mathrm{37}}^{K^{\prime }}$ with temperature (Theroux et al., 2010; Toney et al., 2010). Moreover, the influence of freshwater input on salinity may alter the main alkenone-producing communities in coastal regions (Fujine et al., 2006; Harada et al., 2008) and thus affect the reliability of SST estimates based on the open ocean calibration specifically adapted for Group III alkenone producers (e.g. Emiliania huxleyi). Indeed, strong temperature fluctuations of 10 ^∘C in a Holocene $U_{\mathrm{37}}^{K^{\prime }}$ record from the Sea of Okhotsk were linked to periods with low sea surface salinity, which were in turn correlated to high $U_{\mathrm{37}}^{K^{\prime }}$-derived SSTs (Harada et al., 2008). In contrast, a recent study showed that alkenone producers in particulate organic matter (POM) in a coastal bay in Rhode Island were unaffected by a lower salinity, further illustrated by the excellent match of the 300-year $U_{\mathrm{37}}^{K^{\prime }}$ SST record with instrumental temperature records, despite the proximity to the river (Salacup et al., 2019). Further, Blanz et al. (2005) showed that sediment samples from a salinity transect covering the Baltic Sea to the North Sea showed no relationship between $U_{\mathrm{37}}^{K^{\prime }}$ and SST in the Baltic Proper. Only in the transition zone at Skagerrak, the SSTs were within 1 ^∘C of the global calibration of Müller et al. (1998). Although the high variability in the $U_{\mathrm{37}}^{K^{\prime }}$ SST record and the higher-than-expected reconstructed temperatures at Hank fit with a freshwater input as observed in the Sea of Okhotsk, low BIT index values and T ∕ M ratios in the Hank borehole (Fig. 4) suggest that the organic matter has a primarily marine origin. In addition, the absence of the C_37:4 alkenone in the Hank sediments, a biomarker tentatively linked with coastal or freshwater haptophytes (Cacho et al., 1999), suggests that the $U_{\mathrm{37}}^{K^{\prime }}$ should mostly represent SSTs, although studies from the Baltic Sea indicate that the relative contribution of C_37:4 alkenones only increases at salinities lower than 8 psu (Schulz et al., 2000; Kaiser et al., 2017). Thus, the moderate relation between the percent of C₃₂ diol and $U_{\mathrm{37}}^{K^{\prime }}$ derived SST for the tuned interval (n=26; R²=0.32), suggests that freshwater input may at times have influenced the $U_{\mathrm{37}}^{K^{\prime }}$ SSTs.

Alternatively, the higher $U_{\mathrm{37}}^{K^{\prime }}$ SSTs can be a result of increased production in the spring or summer (Chapman et al., 1996; Rodrigo-Gámiz et al., 2014). Indeed, summer temperatures in the Oosterhout Formation (Ouwerkerk, Netherlands) and contemporaneous Lillo Formation in Belgium (Valentine et al., 2011) recorded from benthic bivalves range from 14.9 to 20.4 ^∘C, which is similar to the range of $U_{\mathrm{37}}^{K^{\prime }}$ SSTs in Fig. 2e. This would mean that summer SSTs were high and very variable during the Pliocene. Although quite variable in the earlier (∼3250–3150 ka) part of the record, $U_{\mathrm{37}}^{K^{\prime }}$ SSTs warmed by approximately 10 ^∘C over the latter part of the tuned interval from 3150 to 3050 ka (Fig. 6d). Reconstructed modern SSTs from surface sediments in Skagerrak region near the opening to the Baltic Sea range from 10 to 12 ^∘C, slightly higher than observed annual SSTs and resembling those of May to June more (Blanz et al., 2005). Thus, there is evidence in the modern North Sea adjacent area for the $U_{\mathrm{37}}^{K^{\prime }}$ recording summer temperatures (coinciding with haptophyte blooms).

LDI SSTs are at first 2 ^∘C cooler than the TEX₈₆ record and then increase toward the same absolute SSTs as in the $U_{\mathrm{37}}^{K^{\prime }}$ record (Figs. 2, 6d). Large discrepancies of 9 ^∘C between TEX₈₆ and LDI-derived SSTs have been observed in the Quaternary of south-eastern Australia (Lopes dos Santos et al., 2013), which the authors attributed to seasonal production of isoGDGTs in the cooler months and long-chain diols in the warmer months. In late Pliocene sediments from the central Mediterranean, LDI SST estimates were slightly lower than $U_{\mathrm{37}}^{K^{\prime }}$ SSTs; however this was within the error range of the proxies (Plancq et al., 2015). Due to the recent advent of the LDI proxy and the scarcity of other multi-proxy studies (De Bar et al., 2018; Lattaud et al., 2018) comparing the LDI to $U_{\mathrm{37}}^{K^{\prime }}$ and TEX₈₆ SSTs in the same sediment samples, further discussion on this topic is limited.

4.4.2 mPWP climate

The mPWP is almost entirely covered by the oxygen isotope age-tuned interval of the Hank record, which starts after the hiatus that marks the M2 event. This mPWP interval is correlated with the Poederlee Formation and Oorderen Sands Member (of the Lillo Formation) located in Belgium (De Schepper et al., 2009a; Louwye and De Schepper, 2010). The average TEX₈₆ (10 ^∘C) and $U_{\mathrm{37}}^{K^{\prime }}$ (16 ^∘C) reconstructed SSTs over the mPWP interval show good agreement with the PRISM3 model reconstructions for February (10.4 ^∘C) and August (16.7 ^∘C; Dowsett et al., 2009). A common feature of the $U_{\mathrm{37}}^{K^{\prime }}$ and LDI SST records is the gradual warming between ∼3150 and 3050 ka (Fig. 6d), seen most clearly in the LDI record. Before the SST warming from 3150 to 3050 ka, the percent of C₃₂ diol decreases slowly (Fig. 6e), indicating a decrease in freshwater discharge and/or an increased distance to the coast. The low T ∕ M ratios and the presence of a clay layer from 292 to 271 m in Fig. 4c (corresponding to 3155–3053 ka; Fig. 6) at this time further indicate increased marine influence, likely as a result of sea level rise. Differences in the absolute degree of warming recorded by the $U_{\mathrm{37}}^{K^{\prime }}$ and LDI SST proxies could be attributed to the different seasons in which they are produced, as well as by lateral transport of certain biomarkers (Benthien and Müller, 2000; Ohkouchi et al., 2002). For example, the change in currents in the North Sea after the M2 event, bringing in warmer waters from the North Atlantic, may have brought alkenones and/or diols with a warmer signature to the SNSB. This would then further contribute to the high SSTs reflected by the $U_{\mathrm{37}}^{K^{\prime }}$ and LDI proxies compared to those recorded by the TEX${}_{\mathrm{87}}^H$ (Fig. 6).

Notably, where the TEX${}_{\mathrm{86}}^H$-derived SST record suggests relatively stable winter temperatures, the LDI and $U_{\mathrm{37}}^{K^{\prime }}$ SST records reflect highly variable SSTs during the mPWP (Fig. 6d). Such high variability is also seen in all other currently available $U_{\mathrm{37}}^{K^{\prime }}$ SST records from the North Atlantic (8 ^∘C, Lawrence et al., 2009; 6 ^∘C, Naafs et al., 2010; 14 ^∘C, Clotten et al., 2018; Fig. S3) and has been explained by a change in the strength of the NAC (Lawrence et al., 2009; Naafs et al., 2010) and orbital forcing (Lawrence et al., 2009). In addition, high variability in the $U_{\mathrm{37}}^{K^{\prime }}$ record from the Iceland Sea was linked to the frequent occurrence of spring sea ice cover and ice-free summers linked to freshwater input (Clotten et al., 2018). At Hank, the influence of changes in the direction and strength of the NAC on $U_{\mathrm{37}}^{K^{\prime }}$ SSTs cannot be unambiguously identified, despite the reoccurrence of Operculodinium centrocarpum during the mPWP. Most of the variation may instead be explained by the varying influence of freshwater input from the proto-Rhine–Meuse River system in combination with the relatively shallow coastal location of Hank, which makes it sensitive to fluctuations in temperature. The influence of orbital forcing on pacing the variation in the NAC (Naafs et al., 2010) and possibly the environmental conditions in the SNSB require further analysis, which is currently not possible in the Hank borehole due to the short length of the tuned interval.

In contrast to the variable marine climate, the terrestrial climate proxies indicate that climate of land was fairly stable. The presence of Taxodioidae-type pollen (Taxodium, Glyptostrobus) throughout most of the mPWP (Fig. 6c) indicates that land temperatures were generally not low enough for prolonged winter frosts. Minimum Taxodioidae-type pollen abundance of 10 % has been associated with a mean temperature of the coldest month of > 5 ^∘C (Fauquette et al., 1998). Both the MAT record (Dearing Crampton-Flood et al., 2018) and the Taxodioidae-type pollen covary throughout the whole record (Fig. 4d), and the absolute MATs and the increased proportion of Taxodioidae-type pollen in the mPWP interval (Fig. 6) support the presumed relatively stable climate conditions on land (Draut et al., 2003; Lisiecki and Raymo, 2005).

Importantly, the new chronology for the Hank sediments provides an opportunity to correlate the stratigraphy concept of the local (Netherlands) qualitative Pliocene–Pleistocene Taxodioidae-type temperature curves proposed by Zagwijn (1960, 1992). It should be noted that the original terrestrial Pliocene stages as summarized by Zagwijn (1992) have not yet been dated independently, and in the type area of the south-east Netherlands, they likely represent much smaller intervals of time compared to the Hank sequence. Zagwijn et al. (1992) inferred mean July temperatures between 15 and 20 ^∘C for the Reuverian, which were placed approximately between 3.1 and 2.5 Ma, with short-lived cool pulses down to ∼12 ^∘C that can also be recognized in the brGDGT MAT record (Fig. 6f). Maximum Taxodioidae abundance and mean July temperatures in excess of 20 ^∘C were reconstructed for the Brunssumian, placed approximately between 3.4 and 3.1 Ma. These reconstructed summer temperatures compare broadly to the SSTs reconstructed using the (presumably) partially summer-biased $U_{\mathrm{37}}^{K^{\prime }}$ proxy, which range between ∼10 and 25 ^∘C in the tuned interval (Fig. 6d). Further correlation and dating studies in proximal boreholes and marine sediment sequences will aid in deciphering the Reuverian A–C substages and assigning absolute ages to the zonation of Zagwijn et al. (1992).

5 Conclusions

The age framework for the mid-to-late Pliocene of the southern North Sea basin (SNSB) constructed here reveals that the M2 glacial event is represented as a hiatus, confirming interpretations at proximal sites in Belgium and the English North Sea coast. Sea surface temperatures were variable, which may be caused by the sensitivity of the shallow Pliocene North Sea to climate change and the influence of freshwater input on lipid biomarker SST proxies. Nevertheless, the variability in SSTs matches that in all other currently available SST records from the North Atlantic, indicating that the marine realm was highly dynamic during the mPWP, probably as a result of shifting currents caused by a reorganization/diversion of the North Atlantic Current. Our terrestrial multi-proxy climate records show a highly consistent signal between lipid biomarker temperatures and pollen assemblages, which show stable terrestrial temperatures of 10–12 ^∘C and the continued presence of warm-adapted tree species during the end of the mPWP. Importantly, the chronology presented here allows placing earlier terrestrial temperature reconstructions for Pliocene NW Europe (Zagwijn et al., 1992) in time. This indicates that the Reuverian Stage concept, characterized by abundant Taxodioidae and Sciadopitys and rare Sequoia, is dated to ∼3.2–2.8 Ma. Further high-resolution analysis will attempt to resolve and date the Reuverian A–C substages in this marine setting.