Are landscapes buffered to high-frequency climate change? A comparison of sediment fluxes and depositional volumes in the Corinth Rift, central Greece, over the past 130 k.y.

Sediment supply from source regions is a fundamental control on the stratigraphic record (Gawthorpe and Leeder, 2000; Allen, 2008a; Blum and Hattier-Womack, 2009; Armitage et al., 2011; Covault et al., 2011; Allen et al., 2013; Hampson et al., 2014; Michael et al., 2014; Romans et al., 2016). In areas of active extension, landscape and rift evolution models have all shown the importance of tectonics in driving the volume, grain size, and characteristics of sediment supply to depocenters (e.g., Humphrey and Heller, 1995; Allen and Densmore, 2000; Gawthorpe and Leeder, 2000; Densmore et al., 2003; Cowie et al., 2006; Backert et al., 2010). However, the importance of climate in controlling sediment fluxes from catchments to depositional basins over a range of timescales remains contentious (e.g., Collier et al., 2000; Jerolmack and Paola, 2010; Simpson and Castelltort, 2012; Armitage et al., 2013; Godard et al., 2013; Foreman, 2014; Hidy et al., 2014). Field-based and seismic stratigraphic studies have demonstrated that sediment fluxes can change from glacials to interglacials (Collier et al., 2000; Fuller et al., 2009; Sømme et al., 2011; Hidy et al., 2014); however, the magnitudes, trends, and causes are still debated. Moreover, some authors have argued that sediment routing systems may be buffered to high-frequency climate variation (e.g., Castelltort and Van Den Driessche, 2003; Jerolmack and Paola, 2010). Additionally, disagreement between numerical models with different setups and assumptions about whether landscapes are insensitive (Jerolmack and Paola, 2010; Armitage et al., 2013) or sensitive (Simpson and Castelltort, 2012; Godard et al., 2013) to glacial-interglacial climate change highlights the need for model validation by further field studies.

In this study, we address these fundamental issues by examining the effect of high-frequency climate change (10⁴–10⁵ yrs) on a complete source-to-sink system: the active Gulf of Corinth rift, central Greece, which is a tectonically and climatically well-understood study area. This sedimentary basin is underfilled, meaning any changes in volumes are due to sediment supply and not accommodation changes. In addition, the basin is closed at glacial lowstands and effectively closed at highstands (Perissoratis et al., 2000), enabling mass balance comparisons. Therefore, this setting allows us to make a simpler comparison of source to sink dynamics than many other settings. A rift naturally produces a series of drainages that wholly drain into the basin, unlike many other types of basins, which are characterized by more complex external and internal drainages. Rift drainages are also relatively small, unless a regional-scale river has been captured by the rift system at one of the ends, which is not the case here. In this contribution, we: (1) predict Holocene sediment fluxes to the basin using the BQART methodology (Syvitski and Milliman, 2007); (2) evaluate the extent to which these fluxes match with depositional volumes calculated from seismic reflection data; and (3) use our seismic constraints on sediment volumes and rates to investigate the response of the Gulf of Corinth system to climate change since 130 ka. The results provide new insights into the sensitivity and response of sediment routing systems on a regional scale to climate over a full glacial-interglacial cycle.

Sediment routing systems are the link between the erosional source and the depositional sink. They allow source and sink to be connected spatially and temporally, enabling mass balance comparisons to be made where sediment budgets can be closed (Allen, 2008a; Covault et al., 2011; Sømme et al., 2011; Allen et al., 2013). Multiple studies have improved our ability to predict and understand the effects and time scales of tectonics on sediment routing systems, suggesting that tectonics can control the magnitude of sediment supply (Sinclair et al., 2005; Cowie et al., 2006; Armitage et al., 2011; Michael et al., 2013), the locus of sediment export (Gawthorpe and Leeder, 2000; Zelilidis, 2000; Cowie et al., 2006; Duffy et al., 2014), and the grain size (Whittaker et al., 2007, 2010). However, there is no consensus on the extent to which erosional-depositional systems are sensitive or buffered to higher-frequency climate change over time periods of 10⁴ to 10⁵ yr (Allen, 2008b).

Numerical modeling studies have investigated the effect of climate on sediment supply from catchments, with differing outcomes, predominantly due to differences in modeling approaches (e.g., Jerolmack and Paola, 2010; Simpson and Castelltort, 2012; Armitage et al., 2013; Godard et al., 2013). The catchment fan model of Armitage et al. (2013) demonstrated that sediment fluxes initially respond to a climate forcing (e.g., the onset of Northern Hemisphere glaciation), but if rainfall continues to oscillate over time scales >100 k.y., the sediment flux response decays rapidly and is therefore buffered. They concluded that initial step changes in climate might be recorded in stratigraphy, but ultimately large sediment routing systems (300–1000 km length) are insensitive to high-frequency precipitation variations (cf. Castelltort and Van Den Driessche, 2003; Allen, 2008b). However, Simpson and Castelltort (2012) argued that catchment systems are sensitive to high-frequency changes in water discharge, which can be related to precipitation changes. Additionally, their model suggested that time-averaged sediment fluxes misrepresent the actual variability of sediment fluxes. The model of Godard et al. (2013) also found similar results to Simpson and Castelltort (2012), where sediment fluxes were demonstrably sensitive to high-frequency changes in precipitation. Finally, a number of modeling studies have identified autogenic behavior in sediment yields, meaning fluxes vary even when there is no change in the system boundary conditions (Jerolmack and Paola, 2010; Van De Wiel and Coulthard, 2010). Such findings therefore suggest that it may not actually be possible to link specific allocyclic forcings to sediment fluxes or volumes; this has been referred to as “signal shredding.”

Clearly, these models do not replicate the full behavior of real-world systems. Consequently, increasing numbers of studies have attempted to address landscape sensitivity to climatic perturbations from field observations. Foreman (2014) re-evaluated the Clark Fork Sheet sandstone in NW Wyoming using δ¹³C records and concluded deposition was the direct result of climate change associated with the Paleocene-Eocene Thermal Maximum, agreeing with both Simpson and Castelltort (2012) and Armitage et al. (2013) models that a step change in climate could lead to an increase in sediment fluxes.

Collier et al. (2000) derived sediment fluxes from sediment volumes in the Alkyonides Gulf, which is part of the Corinth rift (Fig. 1), for the last 130 k.y. and concluded that sediment fluxes based on seismic reflection–derived volumes were ∼60% higher during glacial periods than during interglacials, suggesting that landscapes are sensitive to climate perturbations at a time scale of ∼100 k.y. For insightful observations to be made, studies need to fully investigate the link between catchment erosion and basin deposition, thus investigating the full source-to-sink system. For example, Sømme et al. (2011) investigated the role of both climate and tectonics on the Golo source-to-sink system, Corsica, since 80 ka. They compared BQART-modeled sediment fluxes (Syvitski and Milliman, 2007) with depositional volumes on the shelf. They showed that sediment fluxes were highly sensitive to climate, but the Golo system was buffered on longer tectonic time scales. Similarly, for a series of small catchments in California, Covault et al. (2011) compared cosmogenic radionuclide–derived catchment-averaged denudation rates with seismically derived sedimentation rates offshore and showed that, during falling and lowstand sea level (40–13 ka), sediment budgets balanced well with denudation rates, whereas during rising and high sea levels, sediment budgets were greater than denudation rates.

Strath terraces and paleo-erosion-rate estimates along the Eel River, California, suggest that high erosion rates are correlated to the 70–18 ka glacial period (Fuller et al., 2009), with higher inferred precipitation rates for that region (Adam and West, 1983). However, spatially averaged catchment-wide cosmogenic denudation rates from river catchments in the southern United States implied that interglacial denudation rates were higher than glacial rates (Hidy et al., 2014); these results agreed well with suspended sediment fluxes derived from the temperature-driven BQART model. These results suggest that higher denudation rates were driven predominantly by changes in temperature, a fundamentally different climate regime to that considered by Fuller et al. (2009), and this in turn implies that the lower glacial sediment fluxes in lowland catchments would have been counterbalanced by higher sediment fluxes in mountainous regions (cf. Herman et al., 2013).

Studies of natural systems have therefore provided evidence that particular catchment systems can either respond to or be buffered to short-term climate fluctuations, and that sediment fluxes can be enhanced in either glacial (e.g., Collier et al., 2000; Covault et al., 2011) or interglacial conditions (e.g., Hidy et al., 2014). However, at a regional scale, when integrating many catchments of different sizes, the extent to which sediment routing systems are sensitive to climate fluctuations over time scales of 10⁴ to 10⁵ yrs still remains unclear.

We therefore test these ideas by quantifying sediment accumulation rates and sediment flux estimates for both glacial and interglacial periods in an area where sediment budgets can be closed. We evaluate whether the system is responding to short-term climate change (10⁴ to 10⁵ yrs), and we constrain the magnitude of the response.

The present-day Gulf of Corinth is ∼120 km long (E-W) and ∼30 km (N-S) wide at its widest (Fig. 1). It is one of the fastest extending rifts in the world, with modern-day global positioning system (GPS) extension rates of 5 mm/yr to 15 mm/yr from east to west (e.g., Roberts and Jackson, 1991; Le Pichon et al., 1995; Clarke et al., 1997; Briole et al., 2000; McClusky et al., 2000). Extension in the region initiated in the early Pliocene, at ca. 5m Ma (Ori, 1989; Ford et al., 2013). The basin is asymmetric, with large N-dipping faults on the southern margin creating large depocenters in this region of the gulf for the last 240 k.y. (Nixon et al., 2016). Uplift estimates vary on the south coast, but rates range from 1 mm/yr to 2 mm/yr for Holocene-constrained rates near Eliki (Fig. 1; see reviews of uplift rates in Bell, 2008; Bell et al., 2009). The current drainage network was established at ca. 400–600 ka, when the rift developed its present-day geometry (Ford et al., 2016).

Within the gulf, bathymetry ranges from 860 m to ∼60 m below sea level at the westernmost end (Fig. 1; Nixon et al., 2016), due to the presence of the Rion-Antirion Sill (Moretti et al., 2003). This sill is interpreted to have been tectonically stable since the late Pleistocene, and thus periodically transforming the Gulf of Corinth into a lake during glacial lowstands (Perissoratis et al., 2000). The presence of an open basin and lake over the last 130 k.y. means that the basin has been sediment starved, and so any changes in volume are not due to a lack of accommodation space. Marine-lacustrine periodic transformation is confirmed from long piston cores collected offshore (Moretti et al., 2004) and can be correlated to seismic reflection facies variations in the basin and to marine lowstand clinoforms (Lykousis et al., 1998; Collier et al., 2000; Perissoratis et al., 2000; McNeill et al., 2005; Lykousis et al., 2007; Bell et al., 2008). Marine stratigraphy has high-amplitude seismic facies, and lacustrine stratigraphy has low-amplitude seismic facies (cf. Lykousis et al., 1998; Collier et al., 2000; Perissoratis et al., 2000;, McNeill et al., 2005; Bell et al., 2008). These seismic stratigraphic features have been used to develop an age model for the offshore Corinth rift stratigraphy (Nixon et al., 2016).

We identified four broad climatic periods over the last 130 k.y. for the Eastern Mediterranean: (1) the current interglacial (Holocene), 0–12 ka; (2) the pleniglacial, 12–70 ka, encompassing the Last Glacial Maximum (LGM) at 26.5–19 ka (Clark et al., 2009); (3) interstadials and stadials 1 and 2, 70–115 ka; (4) and the last interglacial, 115–130 ka.

Currently, the Gulf of Corinth experiences a warm temperate climate with hot dry summers (Kottek et al., 2006). Holocene mean sea-surface temperature (SST), based on analysis of scleractinian coral Cladocora caespitosa, is ∼16 °C (Royle et al., 2015). Weather station records for 1950–2000, such as the WorldClim precipitation and temperature data sets (see Hijmans et al., 2005 for full study), suggest a mean temperature range of 11 °C to 17.4 °C, with a yearly mean of 14.5 °C for catchments in this study. Mean annual precipitation for catchments is ∼740 mm/yr.

Long-term pollen records that can be used for reconstructing climate covering a complete interglacial-glacial cycle and the Holocene are not available for the Gulf of Corinth specifically, because no deep cores have been collected. However, elsewhere in Greece, there are records that do span the last interglacial-glacial cycle (e.g., Ioannina, Kopais, Tenaghi—Tzedakis, 1994, 1999; Tzedakis et al., 2003, 2006; and in the Aegean—Triantaphyllou et al., 2009; Kouli et al., 2012). The pollen record at Ioannina, situated in the Pindos Mountains in eastern Greece (see inset in Fig. 1), is topographically and geomorphologically most similar to the Gulf of Corinth.

Data from Ioannina suggest that from the last interglacial into the glacial, there was a decrease in the distribution of temperate forests and an increase in dry shrubland and steppe biomes (Tzedakis, 1994; Fletcher et al., 2010), such as Artemisia (herbaceous plants and shrubs) and chenopods (herbs, shrubs, or small trees; Tzedakis et al., 2002). This increase suggests an increase in aridity during the last glacial (Tzedakis, 1999). However, the presence of arboreal pollen throughout the last 130 k.y. suggests that there was sufficient precipitation and warm winter temperatures to sustain populations of trees (Tzedakis et al., 2002).

Three global circulation models (GCMs; Table 1): Community Climate System Model 4 (CCSM4), Model for Interdisciplinary Research on Climate Earth System Model (MIROC-ESM) and Max-Planck-Institute Earth System Model (MPI-ESM-P) describe both the precipitation and temperature at the LGM. They suggest that mean annual temperatures for the Corinth Rift were slightly warm at the LGM (CCSM4: 9 °C; MIROC-ESM: 8.8 °C; and MPI-ESM-P: 9.6 °C), ∼5 °C lower than the modern WorldClim data (14.5 °C). This is in accordance with the continued presence of Quercus (oak), which requires sufficiently warm winter temperatures and moisture availability (Brewer et al., 2002). The Corinth rift experienced minor glaciation in a few of the catchments (Hughes and Woodward, 2017). SST reconstructions using planktonic foraminifera from Tyrrhenian Sea cores suggest a decrease to ∼10 °C for marine isotope stage 4 (ca. 70–55 ka; Kallel et al., 2000). The GCMs indicate a similar range of precipitation at the LGM as in modern times (CCSM4: 828 mm/yr; MIROC-ESM: 716 mm/yr; and MPI-ESM-P: 977 mm/yr; WorldClim: 925 mm/yr; Table 1).

Pollen records show a series of cycles of open vegetation–steppe (cooler/drier) and sustained tree expansion (warmer/wetter) periods, suggesting that the climate was variable between 115 and 70 ka (Tzedakis, 1994, 1999). The overall percentage of Carpinus betulus (hornbeam) suggests that the climate was warm enough to sustain their population, but the overall decrease in pollen abundance from the last interglacial suggests a cooler climate (Tzedakis et al., 2002). SST reconstructions using planktonic foraminifera from Tyrrhenian Sea cores suggest a relatively warm SST, with a mean for the time period of 14 °C, similar to modern temperatures (WorldClim; Table 1). The proportion of Fagus (beech; Tzedakis et al., 2002) indicates moisture availability persisted long enough for Fagus to become dominant (Roucoux et al., 2011). Overall, approximately one third of this period of time (111–104 ka and 88–83 ka) had colder climatic conditions, and two thirds of the time (104–88 ka and 83–68 ka) had conditions closer to the Holocene and the last interglacial (Tzedakis et al., 2002).

The last interglacial in southern Europe had a duration of ∼15 k.y. During the last interglacial, there was a similar proportion of arboreal pollen to the Holocene (Tzedakis et al., 2002), suggesting similar climates. Mean Tyrrhenian Sea SSTs indicate that the last interglacial climate was slightly warmer (1 °C) than the Holocene (Kallel et al., 2000).

Analysis of a 30-m-spatial-resolution digital elevation model (DEM) from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) was undertaken with the ArcGIS 10.1 hydrological toolbox to derive a drainage network for the Gulf of Corinth. A threshold of 300 pixels (0.27 km²) was used to establish the channel network, and stream order was defined using the Strahler method (Strahler, 1957). We obtained 73 catchments larger than 5 km² draining into the gulf (Fig. 1) and excluded those with an area less than 5 km². These catchments were compared carefully to their geometries on published maps and aerial imagery. A manual adjustment was subsequently made to five catchment boundaries, to bring them into line with our observations, correcting for any artifacts due to the DEM spatial resolution. The sediment flux of each catchment was estimated using the BQART model (Syvitski and Milliman, 2007). We assumed that the drainage network geometry has not changed markedly over the last 130 k.y.

To convert sediment loads from million tonnes per year to m³/yr of sediment, a surface porosity (φ) of 54.8% was used, based on the porosity-depth model of Nixon et al. (2016) for sediment in the gulf, which was derived from the porosity data of Goldhammer (1997). This porosity is consistent with the MD01-2479 long-piston core (Moretti et al., 2004). We used a grain density of 2550 kg/m³, based on an initial bulk density of 1700 kg/m³. Subsequently, we quote the amount of sediment released from the onshore catchments in units of volume per unit time. Consequently, in this study, we defined the term “sediment flux,” a term that is widely used in the source-to-sink community, as the sediment load when in units of m³/yr.

We used the model to estimate Q_s values for the Holocene, and then for the time periods identified in the climate review. Catchment maximum relief and drainage areas were derived from the DEM. For the 12–70 ka period, we added 60 m to the maximum relief, because the Gulf of Corinth is lake during glacial lowstands, rimmed by the Rion Sill (at 60 m below current sea level). We did not change catchment areas in the model between the highstand and lowstand (–60 m) conditions, because (1) the southern margin of the gulf has steep bathymetry offshore (Fig. 1; cf. Nixon et al., 2016), creating little change in plan-view catchment area; (2) the northern margin has shallower bathymetry, and catchments may have been larger by several square kilometers; however, if we did increase every catchment’s area by 5 or 10 km², then it would lead on average to small increases of 3%–6% in modeled sediment fluxes, as catchment area in the BQART equation is to the power 0.5 (unlike relief); and (3) bathymetry and seismic lines do not extend completely to the coast in any case, so we cannot conclusively know the shape of the catchments at the lowstand. We set B equal to 1 because lithological effects on bedrock erodibility are currently not fully constrained (Sklar and Dietrich, 1998) and to avoid unknown variables in our modeling; furthermore, the Gulf of Corinth lacked major glaciation during the glacial period (Hughes and Woodward, 2017). Water discharge estimates for individual catchments were derived by multiplying each catchment area by the mean precipitation estimates for each catchment, as described below. We used this method to estimate water discharge as we can robustly apply it to the last 130 k.y., as we do not have information on runoff:precipitation ratios for the last 130 k.y. Modern data for one river in the area suggests that the runoff ratio is ∼60% (Kaleris and Langousis, 2017); this would reduce modeled sediment fluxes by ∼15% if applied to the system as a whole.

Holocene precipitation and temperature values were derived from the high-resolution WorldClim data sets (Hijmans et al., 2005), created from a compilation of weather stations around the world for the period of 1950–2000, and interpolated with a 30 arc second spatial resolution (Hijmans et al., 2005). The mean temperature for the catchments yields a value of 14.5 °C for the Holocene (Table 1). This temperature is ∼1.5 °C lower than Corinth SSTs (Royle et al., 2015), so we investigated the possibility that Holocene temperatures could be ±10% from the WorldClim values. Annual catchment precipitation means derived from WorldClim data sets range between 544 and 880 mm/yr (Table 1); we similarly considered a ±10% variation from these precipitation values in our model.

Three sets of GCMs for the LGM were used from the Coupled Model Intercomparison Project phase 5/Palaeoclimate Modeling Intercomparison Project phase 3 (CMIP5/PMIP3) representing ∼21 k.y. (Taylor et al., 2012). The LGM models are the Community Climate System Model (CCSM4; Gent et al., 2011), MIROC-ESM (Sueyoshi et al., 2013) and MPI-ESM-P (e.g., Jungclaus et al., 2013). We averaged the three LGM GCMs to produce a synthetic climate model, herein named AV-LGM, to produce end-member sediment flux and subsequent volume estimates for the period of 70–12 ka. We also explored the potential of the climate being 200 mm/yr drier based on our climate review. Finally, we used the MPI-ESM-P LGM GCM data set to estimate sediment fluxes at the LGM, because this is the warmest and wettest of the three LGM GCMs.

For the time period of 115–70 ka, we assumed that the WorldClim data set is representative of the interstadials, which account for two thirds of this time period (see climate review), and the AV-LGM model is representative of the remaining stadial climate conditions. We explored three different setups, where WorldClim and the AV-LGM data set were combined with different ratios: (1) 66:34, (2) 50:50, and (3) 75:25, respectively. For the last interglacial (130–115 ka), because it was climatically similar to the Holocene, we reapplied the modern WorldClim data sets. However, we also investigated the possibility of the climate being warmer by adding 1 °C to each pixel in the WorldClim data set.

Data from four two-dimensional (2-D) seismic reflection surveys were used: the RV AEGAEO 1995; the M.V. Vasilios 2003; the M.V. Vasilios 1996; and the R/V Maurice Ewing 2001. The R/V AEGAEO data set was collected by the Hellenic Centre for Marine Research (here called the HCMR survey) in 1995 (Lykousis et al., 1998, 2007). The survey was collected using a variety of sources (3.5 kHz acoustic profile and 1, 5, 10, and 40 in.³ air guns; Lykousis et al., 1998, 2007). Maximum penetration depth of the synrift sediments was ∼0.4 s two-way traveltime (TWTT; ∼1.25 km).

The M.V. Vasilios 2003 data set was collected by the Universities of Southampton, Leeds, and Patras (McNeill et al., 2005; Bell et al., 2008), and it consists of high-resolution multichannel 2-D seismic reflection profiles collected using a 150–2000 Hz sparker source (Fig. 1; McNeill et al., 2005; Bell et al., 2008). The profiles imaged to depths of 1.5 s TWTT below the seafloor (∼0.75 km).

To aid interpretation within the Alkyonides Gulf (Fig. 1), we used 2-D single-channel sparker data from the M.V. Vasilios 1996 survey, which imaged to 0.85 s TWTT (∼430 m; Fig. 1; Collier et al., 2000; Leeder et al., 2002). Furthermore, we integrated interpretations of the 130 ka seismic horizon from Nixon et al. (2016) that was derived from the deeper, but lower-resolution R/V Maurice Ewing 2001 data set (Taylor et al., 2011) as well as the M.V. Vasilios 1996 survey and the HCMR survey.

Five long piston cores (up to 30 m in length) were collected by the RSS Marion Dufresne in the central rift basin (Moretti et al., 2004). Three cores penetrated through the marine succession into the lacustrine unit, and dating of the sediment at the base of the marine succession for the MD01-2480 core yielded a date of 13,210 ± 50 yr B.P. (Moretti et al., 2004), confirming that the transition from lacustrine to marine conditions was broadly synchronous with the onset of the Holocene (ca. 12 ka).

Seismic stratigraphic techniques (such as downlap, toplap, and onlap terminations) were used to map four seismic horizons in TWTT. These horizons were (1) the seabed, (2) Horizon 1, constrained by piston-core data to be 12 ka, (3) Horizon 2, a distinctive high-amplitude, continuous horizon of unconstrained age, and (4) Horizon 3, a horizon that separates high-amplitude from low-amplitude seismic reflections and is widely accepted to be 130 ka based on seismic facies characteristics consistent with a marine-lacustrine transition (Fig. 2; Lykousis et al., 1998; Collier et al., 2000; Bell et al., 2008; Nixon et al., 2016). Horizon 2 and 3 were not mapped outside the Central, Alkyonides, and Lechaion Gulfs (Fig. 1) due to poor imaging at depth and are therefore only constrained in a dashed-white polygon that can be seen in later in Figure 5. These horizons define three units, the gulf-wide Holocene Unit (seabed–Horizon 1), Unit II (Horizon 1–Horizon 2), and Unit III (Horizon 2–Horizon 3; Fig. 2). In order to make a comparison with the Holocene Unit and Units II and III later, we clipped the Holocene Unit to the white-dashed polygon (seen later in Fig. 5) and named this Unit I. The MD01-2477 long-piston core was correlated to the seismic reflection data in order to determine which seismic reflection corresponds to the most recent marine-highstand flooding event (at ca. 12 ka; Fig. 2); elsewhere in the gulf, we used seismic stratigraphic indicators to interpret Horizon 1 (Fig. 2). Units II and III were defined using the seismic facies principles identified by Lykousis et al. (1998), Collier et al. (2000), Leeder et al. (2002), McNeill et al. (2005), and Bell et al. (2008); the seismic characteristics of these units vary markedly and are explored in the results, below.

Two-way-time structure maps were produced from the horizons bounding the three units and depth converted by applying a linear velocity model developed for the offshore Corinth rift by Nixon et al. (2016), where velocity increases at 1.5 km^–1 s^–1 from the seafloor, where the velocity is 1.55 km^–1 (for further details, see Nixon et al., 2016). The model is based on an integration of tomography (Zelt et al., 2004), prestacked depth-migration velocity models (Clément, 2000), semblance plots (Bell et al., 2008), and geophysical core logs (Collier et al., 2000; Moretti et al., 2004).

To generate a Holocene isopach that covered the entire gulf, we merged our isopach covering the gulf east of Selinoutas (Fig. 1) with a depth-converted version of the isopach from Beckers et al. (2015) for the western gulf (where they assumed average velocity of 1600 ms^–1). We decompacted this whole-gulf Holocene isopach using Equation 2. The volumes of the decompacted Holocene unit, Unit II, and Unit III were additionally calculated for each time interval. We used these volumes to derive sediment accumulation rates, as this methodology helped us to remove the well-known temporal bias associated with deriving sediment accumulation rates from one-dimensional data, such as borehole data (Sadler and Jerolmack, 2015).

Using the BQART model, we estimated that the total yearly sediment flux into the Gulf of Corinth is ∼1,770,000 m³/yr, (i.e., 1.77 km³/k.y.) for our preferred model output using the climate and relief constraints outlined in the section on “BQART Sediment Flux Model.” If these sediment fluxes are assumed to be representative of the Holocene, we predict a total volumetric export of 21.3 km³ in the last 12 k.y. (Table 2). These cumulative sediment fluxes to the Gulf of Corinth could be as high as 1.95 km³/k.y. if mean annual catchment temperatures are 10% warmer than the WorldClim data set, or they could be as low as 1.6 km³/k.y. if Holocene temperatures are colder than the WorldClim data set by 10%, giving plausible Holocene volumes in the range of 19.2 km³ to 23.4 km³. Changing the precipitation rate from our preferred model by ±10% gives values that lie within this range as well (Table 2).

Figure 3 reveals the spatial distribution of BQART-derived suspended sediment fluxes exported into the gulf when using the WorldClim data. Sediment flux estimates are highly variable in space and range; they vary from ∼1000 m³/yr to >300,000 m³/yr, with a median of 9582 m³/yr (Fig. 3). The highest sediment input catchment is the large E-W–striking Mornos catchment on the north coast (Fig. 3); previous studies have shown high sediment outputs from this catchment throughout the last 400 k.y. (Ford et al., 2016). However, apart from the Mornos and two catchments draining into the Gulf of Itea (Fig. 3), the rest of the north coast contribution is a minor part of the overall sediment flux budget.

The south coast catchments are predicted to produce the majority of the sediment supplied to the offshore basin. Catchments that drain the footwall of the West and East Eliki fault, the onshore segment of the Derveni fault, and the Lykoporia fault have the highest sediment fluxes, with the largest exports (Fig. 3): ∼135,600 m³/yr (Selinoutas), ∼97,100 m³/yr (Vourikos), ∼60,000 m³/yr (Krathis), and ∼68,800 m³/yr (Sithas). The Alkyonides Gulf has few large catchments and as a result little sediment input. These sediment fluxes can be converted into estimated sediment yields (i.e., denudation rates) that range between ∼0.18 and 0.55 mm/yr by dividing their respective sediment fluxes by catchment area.

A power-law relationship between modeled sediment flux and catchment drainage area is observed, which is explicit in the BQART equation, with an exponent of ∼1 (Fig. 4A). This exponent value means that suspended sediment fluxes increase almost linearly with area, and therefore the gradient of this linear relationship is a prediction of the average denudation rate for catchments draining into the gulf. From this, we derive a denudation rate, equating to a sediment yield, of 0.36 mm/yr for the modern day catchments.

Figure 4B shows (1) the normalized cumulative frequency (i.e., percentage) of the number of catchments with respect to their individual sediment fluxes (squares), ordered from smallest to largest; and (2) cumulative frequency of the total sediment flux for each of these catchments (circles). The cumulative frequency of the number of catchments (squares) shows that sediment fluxes are not evenly distributed; ∼93% of catchments produce less than 70,000 m³/yr (Fig. 4B). The catchment with the median sediment flux exports just over 10,000 m³/yr. The cumulative frequency of sediment fluxes for these catchments shows that ∼5% of them produce 40% of the total annual budget (Fig. 4B). These 5% of catchments are the largest draining into the gulf, e.g., Mornos and Selinoutas (Fig. 3).

The gulf-wide Holocene Unit, Unit II, and Unit III vary markedly in their seismic character and expression (Fig. 2). The gulf-wide Holocene Unit, 12–0 ka (known as Unit I when clipped to the extent of Units II and III, rather than covering the whole gulf), typically has high-amplitude, high-frequency, continuous, parallel uniform reflections (Fig. 2). A primary diagnostic feature of this unit is the draping and onlapping of reflections onto clinoforms within Unit II around the gulf (Bell et al., 2008). The base is core-constrained to the end of the Pleistocene.

Unit II is highly variable in seismic character spatially and with depth. On the shelf margins and in the Alkyonides, Lechaion, and Eratini subbasins, there are well-developed lowstand clinoforms within this unit (Lykousis et al., 1998, 2007; Collier et al., 2000; Perissoratis et al., 2000; Leeder et al., 2002; McNeill et al., 2005; Bell et al., 2008). To aid the interpretation of Unit II in the central Corinth basin, we separated it into two clear subunits with clear seismic facies differences, IIa and IIb (Fig. 2). IIa ranges from high-amplitude, high-frequency, continuous, parallel uniform reflections to lower-frequency reflections with depth (Fig. 2). A thin package of high-amplitude reflections is observed approximately midway by depth (TWTT) through the unit (Fig. 2), separating Subunits IIa and IIb. This package can be seen throughout the gulf and may represent a short marine incursion at ca. 50 ka (see below), given that high-amplitude reflections are typically interpreted to represent marine sediment accumulation (Collier et al., 2000; Bell et al., 2008; Nixon et al., 2016). Below this thin package bisecting Unit II, reflections have a lower amplitude and lower frequency than they had above in Unit IIa (Fig. 2). Unit III is predominantly high-amplitude reflections throughout, although thinner, continuous low-amplitude areas are present (Fig. 2). The base, Horizon 3, shows a clear transition from the high amplitudes of Unit III to low-amplitude reflections below, which has been interpreted by several authors to be the marine-lacustrine transition at 130 ka (Lykousis et al., 1998; Collier et al., 2000; Bell et al., 2008; Nixon et al., 2016).

While the ages of Horizons 1 and 3 are widely accepted, the age of the base of Unit II (Horizon 2) is more poorly constrained. We interpret a band of high-amplitude reflections within Unit II to be ca. 50 ka (Fig. 2). We base this age on the Siddall et al. (2003) sea-level curve, as a brief marine transgression occurred at that time, and the fact that Collier et al. (2000) found further marine facies at the base of their core older than 51 k.y. B.P., from radiocarbon dating. If the Rion-Antirion Sill has remained tectonically stable at a depth of 60 m, the sea-level curve of Siddall et al. (2003) suggests that lacustrine conditions would have initiated by 70 ka at the latest (Fig. 2). There are earlier points on the Siddall et al. (2003) sea-level curve where the sill depth is approached, and the transition to predominantly lake conditions may plausibly have initiated between 100 ka and 115 ka (Fig. 2). We propose either 100 or 115 ka as realistic ages for the base of Unit II (Horizon 2). We therefore proceed with our age model (Horizon 1 = 12 ka, Horizon 2 = 100 or 115 ka, and Horizon 3 = 130 ka), and we consider the uncertainties in the age of Horizon 2 below.

We estimate a total basinwide decompacted volume for the Holocene of ∼35 km³ using the gulf-wide Holocene isopach (Fig. 5). The seismic reflection profiles do not extend all the way to the coastline (they stop within ∼1 km; Fig. 1), and so the Holocene isopach was extrapolated to the coastline. We estimated that this extrapolated volume is ∼6% of the total volume (this does not include the western end of the gulf, which Beckers et al. [2015] covered with their Holocene isopach). The steep southern margin is poorly imaged in seismic reflection data, but minimal storage of suspended sediment export has likely occurred on these slopes. Coarse bed load is, however, deposited on the shelf in large Gilbert-type deltas; this is less important in this study, because the BQART model only predicts suspended sediment flux, and we did not include the Gilbert deltas in our basin volumes. The BQART model prediction of suspended sediment fluxes for the Holocene is ∼21.3 km³, generated using a standard lithological B value of 1; this compares well with the observed volume measured in the basin and is well within the factor of 5 margin of error stated by Syvitski and Milliman (2007; even if we use the 60% runoff ratio for the Holocene and reduce sediment fluxes by 15%, the BQART model still predicts within a factor of 5). In addition, although the majority of the bed load is captured in Gilbert-type deltas that fringe the gulf, some may make it into the central basin. Therefore, the BQART model (which solely estimates suspended load) is expected to provide an underestimate compared to the observed volume in the basin. Nonetheless, based on our data, we conclude that the BQART model provides realistic first-order estimates of sediment fluxes over time scales of 10⁴ yrs.

Our mapped isopachs allow us to compare sediment volumes deposited in the central gulf during 12–0 ka (Holocene isopach), 100/115–12 ka (Unit II), and 130–100/115 ka (Unit III). The calculation of sediment accumulation rates for these time periods enables us to assess whether sediment flux varied between glacials and interglacials. Figure 6 shows the deeper isopachs generated for Units II and III. We then compare the smaller Holocene Unit (Unit I) with Units II and III, which were calculated over the same areas. The decompacted volumes for Units I, II, and III, within the white dashed polygon, are 16, 80, and 34 km³, respectively (Table 3). The volumes are unambiguously not the same and are fully constrained by seismic mapping. However, the velocity model used for depth conversion (cf. Nixon et al., 2016) could affect the results, and thus we investigated how a plausible ±10% in velocity gradient would impact the volumes obtained. We found that the volumes are largely insensitive to the change in velocity model, with Unit III, at 34 km³, having the greatest uncertainty of ±1.6% difference (Table 3). Although Units I, II, and III do not encompass the entire gulf, we assume that the trends in sediment accumulation rates are indicative of the entire gulf for the last 130 k.y. This is because the spatial extent encompassed by the units is an area that has remained, over the last 240 k.y., the major depocenters in the basin (Nixon et al., 2016). Furthermore, we deliberately made the spatial extent of Units I, II, and III the same to be able to compare them on a like-for-like basis. We now explore different scenarios related to the age uncertainties of the base of Unit II (Horizon 2) and produce sediment accumulation rates for these different scenarios.

By dividing the decompacted volumes of each of Units I, II, and III by their estimated age, we derived average sediment accumulation rates, which reflect changing sediment fluxes to the gulf over the last 130 k.y. The factor with the greatest uncertainty for the sediment accumulation rate calculation is the age model used. Consequently, we determined sediment accumulation rates using three different scenarios to explicitly take into account the uncertainties in the age of the base of Unit II. However, the volume and sediment accumulation rate for the Holocene (Unit I) is unchanged throughout the different scenarios because its base (Horizon 1) has been dated by core. The three scenarios are as follows, and rates are shown in Table 4 and Figures 7A and 7B.

No age for the base of Unit II (Horizon 2) is assumed, and the volumes of Units II and III are summed together to avoid any uncertainty in the age of Horizon 2. This provides a volume for the full period 130–12 ka. The sediment accumulation rate for Unit I is then applied to the last interglacial, 130–115 ka, on the basis that it was climatically similar to the Holocene, and this volume is subtracted from the summed volume of Units II and III.

The base of Unit II (Horizon 2) is assumed to be 100 ka, and therefore Unit II represents 100–12 ka, and Unit III represents 130–100 ka. We can take this further and split Unit III into Subunit IIIb (130–115 ka), where we apply the same sediment accumulation rate as the current interglacial to calculate a volume, and subtract this volume from Unit III to calculate the volume and sediment accumulation rate of Subunit IIIa (115–100 ka).

The base of Unit II (Horizon 2) is assumed to be 115 ka, marking the end of the last interglacial period. Therefore, Unit II is 115–12 ka, and Unit III is 130–115 ka.

Scenario 1 has the fewest horizon age assumptions. If, by inference, the last interglacial period had the same sediment accumulation rate as the Holocene, i.e., ∼1.33 km³/k.y., this would result in the glacial period between 115 and 12 ka having a lower time-averaged sediment accumulation rate of 0.91 km³/k.y., a reduction of 32% relative to the Holocene (Fig. 7A; Table 4).

Scenario 2 generates a Unit II, 100–12 ka, sediment accumulation rate of 0.91 km³/k.y., identical to scenario 1. For Unit III (130–100 ka), the rate is ∼1.13 km³/k.y. If we subdivide Unit III into 130–115 ka (Subunit IIIb) and assign the same sediment accumulation rate as the current interglacial (1.33 km³/k.y.), then the remainder of Unit III, Subunit IIIa (115–100 ka), has estimated sediment accumulation rates of 0.94 km³/k.y. (Fig. 7A).

Scenario 3 assumes the base of Unit II is 115 ka, so Unit III correlates to the last interglacial. This results in Unit III (130–115 ka) having a high sediment accumulation rate of 2.27 km³/k.y. when compared to the current interglacial sediment accumulation rate (1.33 km³/k.y.; Fig. 7B). The glacial period (Unit II: 115–12 ka) in this scenario would have a lower rate at 0.78 km³/k.y. than the other scenarios (Fig. 7B; Table 4).

The trends seen over the last 130 k.y. for scenarios 1, 2, and 3 are similar and all generate lower sediment accumulation rates during glacial periods. The reduction is 32%, 32%, and 66% for scenarios 1, 2, and 3, respectively, relative to the last interglacial. We note a model of constant sediment accumulation rate over time is inconsistent with the seismic interpretations and age models presented here. Instead, for all three scenarios, the volumes constrained from the seismic reflection data reveal that there is an observed reduction in sediment accumulation rates from interglacial to glacial conditions, thus proving there is a signal being recorded in stratigraphy and that the signal is not buffered. Our results demonstrate that sediment accumulation rates in the Gulf of Corinth have clearly changed over the glacial-interglacial cycle—and consequently, so have sediment fluxes into the gulf. We interpret this to reflect a climate-driven difference in sediment supply over this period (see Discussion).

Since the simple climatic and topographic variables used in the empirically derived BQART model give a reasonable match to observed sediment volumes in the basin over the Holocene, we therefore used it to estimate sediment fluxes and volumes for the four different climate periods over the last 130 k.y.: 12–0 ka, 70–12 ka, 115–70 ka, and 130–115 ka. These values were then compared to the sediment volumes and sediment accumulation rates constrained from the seismic reflection data. BQART-derived sediment fluxes are higher for the interglacials (12–0 ka and 130–115 ka) than for the glacial periods (Table 2), where sediment fluxes for the two interglacials are ∼1.77 km³/k.y., but could be as high as ∼1.91 km³/k.y. (Fig. 8; Table 2) if temperatures were warmer by 1 °C for 130–115 ka (Kallel et al., 2000). This is the same trend as the sediment accumulation rates calculated from the seismic reflection data for the three scenarios (Fig. 7).

For 70–12 ka, using the AV-LGM model temperature and precipitation constraints, the BQART model predicts a flux of ∼1.09 km³/k.y. (Fig. 8). If we allow the period between 70 and 12 ka to have a climate that may have been drier with precipitation reduced by 200 mm/yr, consistent with some proxy data, we estimate an end-member sediment flux of 1.00 km³/k.y., a decrease of 8% (Fig. 8). We note that this large change in precipitation results in a relatively minor decrease in sediment flux rates. Furthermore, we used solely the MPI-ESM-P LGM GCM to estimate sediment fluxes, because this was the warmest and wettest of the three LGM GCMs. We found that sediment fluxes could have been as high as 1.19 km³/k.y. In addition, we note that glacial catchments could have been slightly larger due to the lower base level at this time; however, we note that sediment accumulation rates were still lower during glacial periods than during interglacials despite this effect.

During the 115–70 ka period, when climate was fluctuating between warm interstadials and cooler stadial conditions, we used different ratios of WorldClim to the AV-LGM model (see section on “BQART Sediment Flux Model” in data sets and methodology). Using setup 1 temperature and precipitation conditions, the BQART model predictions suggest a suspended sediment flux rate of 1.53 km³/k.y. (Table 2; Fig. 8). Setup 2 and 3 produce end-member estimates of 1.41 km³/k.y. and 1.59 km³/k.y., respectively (Table 2). Combining the sediment fluxes volumes for 70–12 and 115–70 ka, we generated an average sediment flux rate of 1.28 km³/k.y. for the overall glacial period of 115–12 ka (Table 2).

Overall, the rates in Table 2 and Figure 8 show that predicted suspended sediment flux rates were highest during interglacials (1.77 km³/k.y.) and lowest during the pleniglacial (70–12 ka, 1.09 km³/k.y.). Additionally, the period between 115 and 70 ka had an intermediate rate compared to the interglacials and glacials (1.53 km³/k.y.). The BQART model results suggest a reduction of 28% in sediment accumulation between the glacial and interglacial periods, similar to the reduction in interglacial to glacial sediment accumulation rates calculated for all three seismic horizon age scenarios (Fig. 6A; Table 4). Consequently, the simple climatic and topographic variables used to derive the BQART model predict a similar flux response to that calculated independently from the analysis of seismic reflection data.

From our seismic reflection data interpretation, we have shown that sediment accumulation rates over the last 130 k.y. have varied (Figs. 7A and 7B; Table 4). Our results indicate that sediment accumulation rates were lower during the last glacial than during interglacial periods for all three sediment accumulation rate scenarios considered. This implies that sediment fluxes from the source catchments have also varied over a glacial-interglacial cycle. Our data require a sediment accumulation rate of ∼2.9 km³/k.y. for the entire gulf-wide Holocene isopach and 1.33 km³/k.y. for Unit I (see Fig. 5). We obtained ∼0.91 km³/k.y. for both periods between 115 and 12 ka and between 100 and 12 ka (scenarios 1 and 2; Fig. 7A; Table 4), which is a reduction of 32% compared to the present and last interglacial; a much larger reduction was obtained for scenario 3. We infer from this that catchments in the Corinth Rift must have a response time of ∼10⁴ yrs for this signal to be observed in the basin.

Our finding that sediment routing systems at a regional scale are sensitive to climate change is not concordant with the model of Armitage et al. (2013), which investigated landscape response to climate oscillations (100 k.y.), because we observed recordable signals in the basin even though our study focused on the most recent of the many Pleistocene glacial cycles. Our results agree with Simpson and Castelltort (2012) and Godard et al. (2013), who also showed that high-frequency climate changes can alter sediment fluxes, and they are consistent with the study of Hidy et al. (2014), which showed that catchments have 30%–35% higher denudation rates during interglacials, a magnitude similar to this study. However, the changes in sediment flux predictions in the model of Simpson and Castelltort (2012) are related to changes in water discharge; similarly, Godard et al. (2013) showed that sediment fluxes are sensitive to high-frequency precipitation changes. Mean annual precipitation (and therefore water discharge) in our study did not change significantly between the WorldClim data set (Hijmans et al., 2005) and our AV-LGM model used to model the four different climate periods over the last 130 k.y. (see “Climate” section).

Several authors have suggested that landscapes may be buffered to high-frequency environmental change when the driving frequency is shorter than the natural response time, or when the autogenic dynamics dominate (e.g., Jerolmack and Paola, 2010; Van De Wiel and Coulthard, 2010). However, the modeling results of Jerolmack and Paola (2010) suggest a correlation to allocyclic forcing can be made when the forcing time scale is longer than the response time of the system or if the environmental perturbation is large. Our sediment volume data suggest that the landscapes around the Corinth Rift likely have response times shorter than 10⁴ yrs. We should also note that the catchments in this study are relatively small (>5 km² to ∼1000 km²), which may suggest that the reason we see a response to high-frequency climate change is because the response time of such small catchments is shorter than the climate signal (cf. Romans et al., 2016). This is in concordance with studies of similar catchment sizes, such as the Golo system, Corsica, which also responds to high-frequency climate change (Sømme et al., 2011), However our study does not specifically address larger river systems (300 to >1000 km), which may be buffered (Castelltort and Van Den Driessche, 2003). We currently do not have the temporal constraints to investigate millennial-scale sediment flux variations. However, constraining rates at this time scale may be feasible in future when the drilling results of IODP (International Ocean Discovery Program) Mission 381 in the Gulf of Corinth are published.

Comparison with the basin-wide Holocene (12–0 ka) sink volume (35 km³) shows that the BQART model effectively predicts the sediment volume for this time period (21.3 km³), within a factor of 1.6, even when we leave the lithological parameter, B, equal to 1. Because of the general agreement between the BQART model and the Holocene seismic volume, our results suggest that simple climate and topographic variables are sufficient to predict Holocene sediment fluxes reliably to first order.

Our modeled sediment fluxes over the last 130 k.y., using simple climate variables, are higher than the seismically derived sediment accumulation rates of Units I, II, and III within the central basin, Alkyonides Gulf, and Lechaion Gulf (Figs. 7A, 7B, and 8); this is because the seismic units are spatially limited. However, the BQART-derived sediment fluxes show a similar trend to the seismically derived sediment accumulation rates. BQART-derived sediment fluxes for the last 130 k.y. clearly show an increase in sediment flux between glacial to interglacial periods. The percentage change in sediment accumulation rates between interglacials to glacials calculated using scenarios 1 and 2 (32%) is similar to the magnitude change derived with the BQART model (28%). The BQART model, therefore, appears to effectively predict the trends, sensitivity, and approximate magnitude of change in sediment fluxes over the last 130 k.y., at least in this setting.

We show that: (1) seismically derived sediment accumulation rates have varied over the 130 k.y., which we interpret to be driven by climate changes over this time period; and (2) simple climate and topographic variables in the BQART model can effectively predict sediment fluxes over a glacial-interglacial cycle. While this study is the first to investigate sediment fluxes into the whole of the Gulf of Corinth, Collier et al. (2000) carried out a more localized seismic interpretation study into changing sedimentation rates in the Alkyonides Gulf (Fig. 1). They found that sedimentation rates had changed between interglacials and glacials, but, significantly, they argued that glacials produced higher sediment fluxes. We have compared our interpreted Horizons 1, 2, and 3 with those of Collier et al. (2000), and, because we picked very similar seismic horizons, any differences between the studies are not due to differences in seismic mapping. Instead, we attribute the difference in our sediment flux calculations to two factors.

First, we note that the velocity model used by Collier et al. (2000) to depth convert the seismic reflection data is slower in the upper ∼150 m of the basin fill and faster in the deeper section when compared to the model we used. This results in Collier et al. (2000) computing a smaller Holocene volume and larger 70–12 ka volume than would be calculated using the velocity model for this study. Collier et al. (2000) noted that their velocity model is not based on any local seismic data or deep well control, whereas the velocity model used here integrated velocity information from piston cores (Moretti et al., 2004), semblance analysis from seismic data processing (Bell et al., 2008), wide-angle velocity studies (Clément, 2000), and seismic tomography data (Zelt et al., 2004), and it is therefore likely to better capture the true velocity structure.

Second, although our horizon interpretations are similar those of Collier et al. (2000), we differ in the assigned age given to the base of Unit II (Horizon 2). Collier et al. (2000) gave this horizon an age of 70 ka; although this is not constrained by core data, it presumes the Rion Sill has remained at a depth of 60 m and uses the Shackleton (1987) sea-level curve. As mentioned previously, water depths would have approached the sill depth, and the first transition to lake conditions may plausibly have initiated at the earlier time of 115–100 ka (Fig. 2).

Although Collier et al. (2000) had no direct age-control on the horizon they estimated at 70 ka, they did have a marine incursion dated older than 51 k.y. B.P., which stratigraphically lies above their 70 ka horizon (or our Horizon 2). Between Horizon 1 and 2, we observed a suite of high-amplitude reflections, which we suggest may be correlated to the marine incursion at ca. 50 ka, based on the Siddall et al. (2003) sea-level curve (Fig. 2). If these high-amplitude reflections are correlated correctly to ca. 50 ka, and we use the Collier et al. (2000) 70 ka age estimate for Horizon 2, we would calculate sediment accumulation rates for our volumes of ∼1.1 km³/k.y. for 50–12 ka (IIa; Fig. 2) and 2 km³/k.y. for 70–50 ka (IIb; Fig. 2). When we compare our 12–0 ka Holocene rate of 1.33 km³/k.y. with the rate calculated for the 70–50 ka period of 2 km³/k.y., it would mean that this interval produced ∼50% higher fluxes than the current interglacial.

The sediment accumulation rates between 70 and 50 ka would actually be the highest sediment fluxes over the entire 130 k.y. period using our calculated volumes. Based on our climate review, we see no reason why the most significant change in sediment accumulation rates over the complete 130 k.y. interglacial-glacial-interglacial cycle should occur between 70 and 50 ka. End-member changes between interglacial and glacial conditions are the most significant climate change over the last 130 k.y.; the 70–50 ka period is a transition between interstadial-stadial conditions to peak glacial conditions in which the amplitude of the climate change and its rate were lower than elsewhere in the glacial cycle. We therefore suggest that it is unlikely that a peak in sediment accumulation rate should correspond to this interval, and we propose an age of 115/100 ka for Horizon 2 to be more likely than the 70 ka age of Collier et al. (2000). In summary, the combination of (1) a different (but probably more realistic) velocity model and (2) a difference in age assigned to Horizon 2 (70 ka vs. 115/100 ka) accounts for the differences between our study and that of Collier et al. (2000). Collier et al. (2000) concluded that glacial periods produced more sediment export, and they inferred that cold wet winters were the driving force behind these higher rates. However, palynology records suggest that glacial periods in Greece are, if anything, slightly drier than interglacials, based on the presence of Artemisia and chenopods (Tzedakis, 1999; Tzedakis et al., 2002). Our sediment flux model indicates higher predicted suspended sediment fluxes during interglacials than glacials, which is consistent with our sediment accumulation rate estimates. When combined with palynology information on climate, we would suggest warmer interglacials drive sediment supply changes. This finding is in direct concordance with the field studies conducted by Fuller et al. (2009) and Hidy et al. (2014), where higher sediment fluxes were correlated to periods of higher precipitation in California during the last glacial and higher temperatures in interglacials, respectively.

Our model results allow us to evaluate the role of temperature and precipitation in driving changes in sediment flux. Mean annual precipitation estimates between the interglacial/modern (WorldClim) and LGM conditions (AV-LGM) used in the BQART model are insufficient to greatly alter sediment fluxes, and so suspended sediment flux change cannot be explained by simple water discharge changes. The reduction in temperature in the BQART model is therefore primarily responsible for the reduction in model-predicted glacial sediment fluxes, although this in itself does not provide a physical mechanism for enhanced sediment fluxes during the warmer interglacials. Next, we consider two alternative hypotheses to explain the observed sediment flux changes.

The comparison between our seismically derived volumes in the gulf and our BQART modeling of sediment fluxes indicates that sediment accumulation in glacial conditions was lower than interglacials. A simple hypothesis would be that glacial conditions result in reduced export due to lower sediment production. Some numerical studies have argued that within Mediterranean catchments, glacial sediment production is lower due to steppe vegetation being dominant (Leeder et al., 1998). However, the same study argued that sediment discharges were actually higher during this period (Leeder et al., 1998), a situation that we do not observe in our basin fill.

The decrease in temperature combined with a similar amount of precipitation predicted for the Corinth rift between the interglacial and glacial periods would typically imply an increase or maintained level of physical weathering during glacial periods (Tucker et al., 2011; D’Arcy et al., 2017). Indeed sediment production is often considered higher in glacial periods due to vegetation changes and periglacial processes such as more effective freeze-thaw weathering (e.g., Matsuoka, 2008) and increased sediment generation on hillslopes. This is supported by the study by Hales and Roering (2007) in the Southern Alps, New Zealand, where they showed that sediment generation rates are higher in present-day alpine conditions above certain high elevations. Furthermore, during the LGM, this altitudinal zone would have been at a considerably lower elevation (Hales and Roering, 2007). Indeed, increased physical weathering in the Mediterranean has been proposed as one reason for the fact that many preserved/exposed bedrock fault scarps postdate the LGM (e.g., Benedetti et al., 2002; Palumbo et al., 2004; Roberts and Michetti, 2004), as fast periglacial erosion rates and resultant high sedimentation rates in the hanging wall are argued to have buried (Benedetti et al., 2002; Roberts and Michetti, 2004) or destroyed (Roberts and Michetti, 2004; Tucker et al., 2011) scarps formed during the peak of the last glacial period.

Consequently, we suggest that the reduction in sediment accumulation rates during the glacial may not be tied simply to changes in sediment generation within catchments, but it may reflect changes in sediment discharge from the catchments to the depocenter. In this case, sediment production may be maintained or elevated, suggesting an out-of-phase situation, where production is high in glacials, but export to neighboring depocenters is lower than during interglacial periods. The recognition that this may be an out-of-phase situation is important within the context of signal preservation. In this case, it would mean that changes in volumes in stratigraphy represent the changing sediment fluxes caused by different climate regimes; the volume changes would not reflect differing rates in primary sediment production (cf. D’Arcy et al., 2017).

An alternative explanation for higher sediment accumulation rates in interglacials is that while mean annual precipitation does not change significantly between glacial and interglacial conditions, the distribution of annual rainfall does. It is already known that storm intensity scales predictably with atmospheric moisture content and temperatures (Trenberth et al., 2003; Berg et al., 2013). We therefore hypothesize that warmer interglacial climates are more conducive to generating higher-intensity storms than colder interglacials (cf. Li and Battisti, 2008; Miller et al., 2010; D’Arcy et al., 2016), generating more threshold-surpassing events that could do more geomorphic work. Storms have been shown to cause higher sediment fluxes directly into the Corinth basin; for example, Heezen et al. (1966) directly correlated breaks in submarine telephone cables that ran along the southern Gulf of Corinth continental slope to submarine landslides directly related to storms. Similarly, Mamassis and Koutsoyiannis (1996) explained that intense storms are associated with summer periods for the modern Greek system. Therefore, although water discharge has relatively little significance in the BQART equation, the variance in water discharge may be taken into account by temperature because increased variability in precipitation increases with higher temperature (cf. D’Arcy et al., 2017). Our hypothesis that variance in water discharge is controlling sediment export is comparable to the numerical modeling results of Godard et al. (2013) and Simpson and Castelltort (2012).

In this model, sediment production in the glacial period, which may have been enhanced relative to today, is initially stored in the catchments surrounding the Gulf of Corinth. Subsequently, in the interglacial, higher discharge events incise, entrain, and export this sediment to the Gulf of Corinth depocenters. Heightened sediment production in glacials agrees with the results of Collier et al. (2000), but we argue that it is stored and “flushed out” in interglacials. Figure 9 synthesizes and visually illustrates this hypothesis.

Studies of Mediterranean alluvial fans have shown that periods of fan aggradation correlate to glacials, where sediment supply to fans increase and periods of dissection relate to interglacials (Harvey et al., 1999; Pope and Wilkinson, 2005), and this fan incision would lead to an increase in sediment being supplied to catchment rivers. Our field observations show that the river systems we studied on the southern margin of the Gulf of Corinth do have incised river-fill terraces that are undated but attributed to the late Pleistocene, and that could plausibly be related to heightened sediment accumulation within the catchment during the last glacial and subsequent incision in the current interglacial (cf. Church and Slaymaker, 1989). We suggest this combination of drivers is the most plausible explanation for the sediment accumulation rate changes observed in the Gulf of Corinth in this study.

Higher-temporal-resolution stratigraphic constraints, currently being generated from IODP Mission 381, will provide further data with which to test these ideas in the future, as sediment accumulation rates will be available over much finer time intervals. In addition, for our hypothesis to be rigorously tested, i.e. that sediment production in glacial periods is heightened, further investigation would be needed on the role of sediment production on hillslopes and their storage capacity over short time periods. Similarly, more work is required on quantifying the amount of incision on river-fill terraces throughout the river systems draining into the gulf to constrain the sediment storage in the catchments at the time of deposition and subsequent incision and sediment fluxes. This would require the dating of river terraces by techniques such as cosmogenic-nuclide analysis or optically stimulated luminescence to constrain their ages. Such data, in combination with terrace mapping, could then be used to better estimate sediment storage within the catchment at the time of river terrace deposition, and then subsequent incision rates could be quantified.

This study compared measured sediment volumes in the Gulf of Corinth basin with predicted sediment fluxes over the last 130 k.y. to investigate how climate controls sediment flux. The principal conclusions are:

• (1) The Holocene sediment volume for the entire Gulf of Corinth, from seismic reflection data interpretation, is 35 km³; the BQART sediment flux model predicts this volume effectively within a factor of 1.6. This result suggests that simple topographic and climatic metrics are sufficient to provide first-order estimates of sediment fluxes over time scales of 10⁴ yrs.

• (2) Sediment accumulation rates in the basin over the last 130 k.y. are ∼30% higher in interglacials than glacials.

• (3) The BQART modeled sediment fluxes predict the same magnitude of decrease from interglacial to glacial periods as observed in the basin sediment volumes, so the simple climatic and tectonic variables used in the model can predict temporal changes in sediment accumulation and sediment fluxes.

• (4) We conclude that landscapes are sensitive to high-frequency (tens of thousands of years) climate change, because the offshore sediment accumulation rates reflect changes between glacial (lower rates) and interglacial (higher rates) periods. These changes must reflect actual changes in sediment fluxes over the last 130 k.y. because the basin is overall underfilled and sediment starved.

• (5) We hypothesize that these changes in sediment accumulation rates are a product of heightened sediment production in glacials due to a combination of lower temperatures (∼5 °C) yet similar mean annual precipitation rates compared with interglacials. This sediment is stored in the catchments until subsequent warmer interglacials, which experience more intense storms. This results in higher-discharge events, which mobilize and export sediment into the gulf.

• (6) The changing volumes and sediment accumulation rates over an interglacial-glacial cycle lead us to the conclusion that catchments are responding to short-term climate signals, which are recorded in the basin. However, we argue that these signals are modified by upstream sediment storage and release, and these effects can complicate the coupling between landscape denudation and observations of basin fill.

This research was funded by the Imperial College Janet Watson bursary, the Natural Environment Research Council (NERC) Science and Solutions for a Changing Planet Doctoral Training Partnership support to Watkins, and a Royal Society grant to Bell and Whittaker. We would like to thank Alex Coleman for his feedback on the decompaction of the seismic volumes, and Patience Cowie for her insightful comments on the discussion. We would also like to thank Demitris Sakellariou and the Hellenic Centre of Marine Research for the use of the R/V AEGAEO seismic reflection data set; Brian Taylor for the R/V Maurice Ewing seismic reflection data set; and Richard Collier, Mike Leeder, Mark Trout, George Ferentinos, Evrivriadis Lyberis, and George Papatheodorou for the M.V. Vasilios seismic reflection data set in the Alkyonides Gulf. Furthermore, we would like to kindly thank the three reviewers, Paul Umhoefer, Brian Romans, and one anonymous reviewer, for their useful comments and suggestions to improve this paper. In addition we would also like to thank Arnaud Beckers for the use of his western Gulf of Corinth Holocene isopach in this study.