Paleo‐denudation rates suggest variations in runoff drove aggradation during last glacial cycle, Crete, Greece

Fluvial aggradation and incision are often linked to Quaternary climate cycles, but it usually remains unclear whether variations in runoff or sediment supply or both drive channel response to climate variability. Here we quantify sediment supply with paleo‐denudation rates and provide geochronological constraints on aggradation and incision from the Sfakia and Elafonisi alluvial‐fan sequences in Crete, Greece. We report seven optically stimulated luminescence and ten radiocarbon ages, eight 10Be and eight 36Cl denudation rates from modern channel and terrace sediments. For five samples, 10Be and 36Cl were measured on the same sample by measuring 10Be on chert and 36Cl on calcite. Results indicate relatively steady denudation rates throughout the past 80 kyr, but the aggradation and incision history indicates a link with climate shifts. At the Elafonisi fan, we identify four periods of aggradation coinciding with Marine Isotope Stages (MIS) 2, 4, 5a/b, and likely 6, and three periods of incision coinciding with MIS 1, 3, and likely 5e. At the Sfakia fan, rapid aggradation occurred during MIS 2 and 4, followed by incision during MIS 1. Nearby climate and vegetation records show that MIS 2, 4, and 6 stadials were characterized by cold and dry climates with sparse vegetation, whereas forest cover and more humid conditions prevailed during MIS 1, 3, and 5. Our data thus suggest that past changes in climate had little effect on landscape‐wide denudation rates but exerted a strong control on the aggradation–incision behaviour of alluvial channels on Crete. During glacial stages, we attribute aggradation to hillslope sediment release promoted by reduced vegetation cover and decreased runoff; conversely, incision occurred during relatively warm and wet stages due to increased runoff. In this landscape, past hydroclimate variations outcompeted changes in sediment supply as the primary driver of alluvial deposition and incision.

samples, 10 Be and 36 Cl were measured on the same sample by measuring 10 Be on chert and 36 Cl on calcite. Results indicate relatively steady denudation rates throughout the past 80 kyr, but the aggradation and incision history indicates a link with climate shifts. At the Elafonisi fan, we identify four periods of aggradation coinciding with Marine Isotope Stages (MIS) 2, 4, 5a/b, and likely 6, and three periods of incision coinciding with MIS 1, 3, and likely 5e. At the Sfakia fan, rapid aggradation occurred during MIS 2 and 4, followed by incision during MIS 1. Nearby climate and vegetation records show that MIS 2, 4, and 6 stadials were characterized by cold and dry climates with sparse vegetation, whereas forest cover and more humid conditions prevailed during MIS 1, 3, and 5. Our data thus suggest that past changes in climate had little effect on landscape-wide denudation rates but exerted a strong control on the aggradation-incision behaviour of alluvial channels on Crete. During glacial stages, we attribute aggradation to hillslope sediment release promoted by reduced vegetation cover and decreased runoff; conversely, incision occurred during relatively warm and wet stages due to increased runoff. In this landscape, past hydroclimate variations outcompeted changes in sediment supply as the primary driver of alluvial deposition and incision.

K E Y W O R D S
alluvial fan, cosmogenic nuclides, fluvial aggradation, incision, paleo-denudation rates, postburial production, sediment supply 1 | INTRODUCTION Quaternary glacial cycles are thought to influence sediment production, transport, and deposition, by changing the transport capacity of streams and the amount of hillslope sediment supplied to them (Bull, 1991). This assumption is based on commonly observed coeval alluvial aggradation and/or incision chronologies that coincide with past climate variations. However, alluvial records respond to changes in hillslope-to-channel sediment delivery and runoff, and it is challenging to disentangle whether a switch to aggradation or incision is driven by changes in sediment input, water discharge, or both . Several Mediterranean studies have reported aggradation of rivers during stadial (cold) stages of the glacial cycle, with one explanation being intensified denudation processes, such as frost cracking, increasing sediment supply (Coltorti & Dramis, 1995;Macklin et al., 2002;Rose et al., 1999;Wegmann & Pazzaglia, 2009). However, paleo-denudation studies on local (e.g. Berger et al., 2008;Shuster et al., 2005;Valla et al., 2012;Vernon et al., 2008) and global scale have produced contradictory results regarding an increase in denudation during glacial climate (Herman et al., 2013;Peizhen et al., 2001;Schildgen et al., 2018;Willenbring & von Blanckenburg, 2010). Linking Quaternary glacial cycles with changes in denudation rates at large scales is complicated by locally variable denudation processes. For example, non-glaciated areas may respond differently to climate change than glaciated regions (Mariotti et al., 2021). Therefore, studying the combined alluvial and denudation response to climate change in controlled settings is an important step towards understanding the controls on aggradation and the impact of climate change on landscapes.
Cosmogenic nuclide concentrations in stream sediment provide a useful tool to quantify hillslope denudation and sediment supply throughout the Quaternary (e.g. Kapannusch et al., 2020;Lenard et al., 2020;Mariotti et al., 2021). Cosmogenic nuclide concentrations from alluvial sediments can be used to calculate paleo-denudation rates, provided concentrations are corrected for nuclide production and decay after sediment deposition (Schaller et al., 2002). These paleo-denudation rates can then be compared with past climate (Mason & Romans, 2018) and aggradation-incision records  to study the geomorphic response to climate change.
However, most paleo-denudation rate studies measure 10 Be in quartz-bearing rocks, creating a knowledge gap regarding the climatic sensitivity of regions with carbonate bedrock.
Several studies have demonstrated that aggradation of alluvial and submarine fans record variations in catchment sediment export in response to Quaternary climate cycles Macklin et al., 2002;Pope & Wilkinson, 2005;Waters et al., 2010;Watkins et al., 2019). However, linking these variations to sediment supply has proven difficult as it remains unclear if denudation in fluvial landscapes increased during a colder climate due to enhanced physical erosion (Hales & Roering, 2007;Marshall et al., 2021). In temperate climates, Marshall et al. (2015) showed a 2.5-fold increase in denudation during stadial intervals using cosmogenic nuclides, whereas other studies do not detect significant variations in denudation over the Pleistocene (Cyr & Granger, 2008;Kapannusch et al., 2020;Schaller et al., 2002). In the southwestern United States, Pleistocene paleodenudation rates show only minor variations across orbitally driven hydroclimate cycles (Mason & Romans, 2018), despite major concurrent variations in catchment sediment flux . As short source-to-sink systems, alluvial fans respond rapidly to changing climatic conditions. When coupled with paleo-denudation rate measurements, alluvial fans are ideal for investigating how climate changes affect hillslope sediment supply to channels, and how such variations in supply are transferred into the sedimentary record (Savi et al., 2014).
To investigate the response of hillslope denudation and stream aggradation-incision behaviour to climate variations in a fluvial landscape, we studied Late Pleistocene deposits of the Sfakia and Elafonisi fan complexes in western Crete, Greece ( Figure 1). We take advantage of geochronologic ages published for the Sfakia fan (Pope et al., 2008(Pope et al., , 2016 and contribute new geomorphic mapping results along with luminescence and radiocarbon dating for the Elafonisi fan system. In addition to studying changes in hillslope denudation, we measured cosmogenic 10 Be and 36 Cl in quartz-bearing and limestone lithologies, respectively, from modern and Late Pleistocene sediment. Together, these data allow us to compare variations in aggradation and incision with concurrent hillslope denudation rates and investigate the sensitivity of the sedimentary record to variations in runoff and sediment supply.  (Reilinger et al., 2006(Reilinger et al., , 2010. Most of Crete is part of a nappe pile stacked during mid-Cenozoic subduction and exhumed during subsequent extensional deformation in the late Cenozoic (Fassoulas et al., 1994;Jolivet & Brun, 2010). The Sfakia fan is a large ($3.2 km 2 ) telescopic fan and the westernmost of seven coalesced fans (Nemec & Postma, 1993). The 24 km 2 catchment sourcing the Sfakia fan consists mostly of platy limestones with interbedded cherts of the Jurassic to Eocene Plattenkalk Formation and limestone of the Jurrasic Trypali unit (Creutzberg et al., 1977). The Elafonisi fan system consists of eight adjacent alluvial fans along the southwestern coast of Crete. We collected samples along the Stomio River, which drains the largest (35 km 2 ) catchment of the Elafonisi area and has several river terraces preserved near its outlet. Additionally, we sampled two $1.5 km 2 catchments with telescopic fans near Livadia and Chrissoskalitissa ( Figure 2 (Shaw et al., 2008) or a sequence of earthquakes  by about 8 m at Elafonisi and 2.9 m at Sfakia (Angelier, 1979;Ott et al., 2021). However, on a Late-Pleistocene timescale, rock uplift rates along the coast of Crete are relatively steady (Gallen et al., 2014;Ott et al., 2019b;Robertson et al., 2019;Strasser et al., 2011).

| BACKGROUND
In contrast to the Elafonisi fans, the Sfakia fan has been studied in detail for its sedimentology and depositional history. Two main units can be distinguished (Nemec & Postma, 1993;Pope et al., 2008Pope et al., , 2016, which we refer to as Quaternary deposit 1 (Qt1) and Qt2. Qt1 is only preserved at the head of the Sfakia fan and is dominated by thick beds of subangular to subrounded, cobble-to-boulder-sized clasts, interpreted as debris-flow deposits (Nemec & Postma, 1993;Pope et al., 2008). Unit Qt2 is inset into Qt1 and consists of bettersorted, pebble-to-cobble streamflow deposits with interspersed debris-flow layers. Both fan units contain layers of fine-grained red silt with angular pebbles and pale tan silt lenses, representing paleosols and small aeolian dune deposits, respectively (Nemec & Postma, 1993). For a more detailed description of the sedimentology of the Sfakia fan, we refer the reader to Nemec and Postma (1993) and Pope et al. (2008). Based on published geochronologic ages, Qt1 was deposited mainly during MIS 5, and after a brief period of incision at the beginning of MIS 4, Qt2 was deposited from MIS 4 to 2 (Pope et al., 2008(Pope et al., , 2016. Qt2 is now deeply incised; based on archaeological evidence, the incision must have started before about 5.5 ka (Pope et al., 2008).
Today, both catchments are dominated by a hot-summer Mediterranean climate (Csa in the Köppen-Geiger climate classification) (Beck et al., 2018). The mean annual precipitation for the Elafonisi and Sfakia catchments is 817 and 656 mm/a, respectively. The lower precipitation rate and reduced water-holding capacity of the limestone terrain dominating the Sfakia catchment likely explains its lower vegetation density compared to the Elafonisi area (Ott, 2020).
Eastern Mediterranean climate during the last glacial cycle was characterized by cold and dry stadials (MIS 2, 4, 6) and comparatively warm and wet interglacials and interstadials (MIS 1, 3, 5) (Cheddadi & Rossignol-Strick, 1995;Langgut et al., 2011). A potential glaciation of the Levka Ori Mountains above the Sfakia catchment remains debated (Hughes & Woodward, 2017). Nevertheless, even if small glaciers existed in the Levka Ori Mountains, the Sfakia catchment presently lacks surficial streams connecting the fan to the highest areas of the Levka Ori. A potential hydrologic connection to the high-altitude glaciated areas may have existed, but due to the small extent of the potentially glaciated area compared to the catchment size (24 km 2 ), we assume no major glacial influence on the Sfakia fan.
Catchment-wide denudation rates on Crete have been determined from catchments draining the PQ bedrock unit in western Crete and carbonate catchments throughout the island (Ott et al., 2019a). 10 Be-derived denudation rates in western Crete are around 0.11 mm/a, and 36 Cl-derived denudation rates from carbonate areas average 0.13 mm/a (Ott et al., 2019a). Ott et al. (2019a) concluded that the carbonate massifs and PQ cored mountains in western Crete denude at similar rates but that infiltration of precipitation into the karst system permits the carbonate mountains to form steeper topography despite experiencing similar uplift rates.

| MATERIALS AND METHODS
We built on the existing geochronological framework established by 39 luminescence and U-Th ages from the Sfakia fan (Pope et al., 2008(Pope et al., , 2016 and collected samples for paleo-denudation rates from dated and rapidly buried layers. For the previously undated Elafonisi fan system, we established the geochronology using luminescence and radiocarbon dating. Samples for paleo-denudation rates were taken from sedimentary layers with age control, and either 36 Cl calcite or 10 Be quartz concentrations were measured, depending on the sediment source.

| Field mapping and digital elevation data
The redness of soil B-horizons in arid landscapes is related to the concentration of secondary iron oxides and has been used extensively for relative age determination of alluvial fan surfaces (Ferrier & Pope, 2012;Harvey et al., 1999). The redness of B-horizons within soils was determined using a Munsell colour chart on the fine-grained matrix between clasts. Mapping of the Quaternary (Qt) units is based on soil redness, geomorphic position (elevation), observed unconformities, and abrupt changes in fan surface slope. At the Sfakia fan, we used a drone for aerial images and built a digital elevation model (DEM) using photogrammetric software (Agisoft). At Elafonisi, a 5 m DEM was provided by the Hellenic Cadastre SA. For mapping marine paleoshorelines, the inner shoreline angle (ISA) was estimated in the field, and its elevation was measured with a laser range finder. ISA refers to the intersection point between the marine abrasion platform and the paleo-sea cliff and approximates the paleo-sea level during paleoshoreline formation (Lajoie, 1986).

| Optically stimulated luminescence
We collected seven samples from fan and river-terrace deposits within the Elafonisi study area to constrain the timing of sediment deposition and burial by optically stimulated luminescence (OSL) measurements (Table 1). Samples were collected in light-proof metal cylinders driven 0.25 m horizontally into sand-dominated sedimentary lenses. In addition, moisture samples were collected in airtight film canisters from the sample cylinder location, and environmental dose rate samples were collected from sediment within an $30 cm radius surrounding the OSL sample.
OSL samples were analysed by single-aliquot regenerative dose (SAR) procedures for dating quartz sand (Murray & Wintle, 2000;Wintle & Murray, 2006) at the Utah State University Luminescence Laboratory, Logan, UT, USA. Samples were processed under safelight conditions using sieving, gravity separation, and acid treatments with HCl and HF to isolate the quartz component of the 90-150 μm grain size range following procedures outlined by Aitken (1998) and described in Rittenour et al. (2003Rittenour et al. ( , 2005. Samples were analysed on Risø OSL/TL DA-20 luminescence readers with blue/green (470 nm, 36 W/m 2 ) stimulation and detection through 7.5-mm UV filters (U-340). Stimulation was conducted at 125 C following 240 C preheats (10 s) for regenerative and natural doses. Dose-response curves were created using five sensitivity-corrected regenerative dose points (including a repeat and zero-dose step). The resultant data were fit with a saturating exponential curve from which the equivalent dose (D E ; reported in Figure S1) is determined. The D E , reported in Table 1, is based on a central age model (Galbraith & Roberts, 2012) from the measurement of at least 20 aliquots of sand mounted on a 2 mmdiameter area of the measurement disks ($200 grains/aliquot). Dose rate measurements were determined by chemical analysis of the U, Th, K, and Rb content using ICP-MS and ICP-AES techniques and conversion factors (Guérin et al., 2011) (Table S4). The contribution of cosmic radiation to the dose rate was calculated using sample depth, elevation, and latitude/longitude following Prescott and Hutton (1994). Sediments were dry when collected, but a value of 3 AE 3 wt% water content was considered representative of their burial history.

| Radiocarbon dating
To constrain the end of the last depositional phase at the Sfakia fan, we dated carbonate rinds (thickness 1-5 mm) from the upper metre of a coastal cliff at the distal end of the Sfakia fan ( Table 2). The dating of carbonate cement has been utilized in similar studies to constrain the timing of cementation periods, thereby providing a minimum depositional age (Geyh & Eitel, 1997;Giresse & Martzluff, 2015). All samples were leached with 0.5 M HCl to remove possible surface contamination, put in 12 ml septum-sealed vials (Labco), flushed with , and converted to graphite targets using automated graphitization equipment (Synal et al., 2007).
Accelerator mass spectrometry (AMS) was carried out at the Laboratory for IonBeam Physics at ETH Zürich using a Mini Carbon Dating System (MICADAS) AMS (Synal et al., 2007). Oxalic acid II NIST standard was used for fractionation correction and standard normalization.
Additionally, we dated a terrestrial gastropod (snail) shell sampled from a silt layer contained within a channel incised into the fan surface at the Chrissoskalitissa site to constrain the timing of fan-surface abandonment and the onset of incision (Table 1, Figure 3G). After sieving bulk sediment and microscopic examination, a single, intact shell was selected. The shell was broken and cleaned in milli-Q water with sonic vibration before removing the outer surfaces with dilute HCl. This sample was analysed at the University of California Irvine Keck Carbon Cycle AMS Facility (UCIAMS), following standard cleaning and preparation techniques. The reported radiocarbon age was corrected for isotopic fractionation according to Stuiver and Polach (1977), and sample preparation backgrounds were subtracted based on measurements of 14 C-free calcite. The sample δ 13 C value was measured relative to standards traceable to Pee Dee Belemnite, using a Thermo Finnigan Delta Plus stable isotope ratio mass spectrometer with gas bench input.

| Cosmogenic nuclide measurements
We analysed 11 cosmogenic nuclide samples to determine (paleo-) denudation rates. For three samples sourced by the quartz-rich PQ unit (Stomio River and Livadia at Elafonisi), we measured only 10 Be on quartz (

| 10 Be measurements
In the northern Elafonisi area, we collected three samples of sand and processed the 0.25-0.71 mm size fraction for 10 Be analysis. We separated quartz from other minerals by magnetic separation and repeated etching with hydrochloric, hydrofluorosilicic, and hydrofluoric acid. A 9 Be carrier was added, and beryllium was extracted using ion chromatography (Bierman et al., 2002). The 10 Be/ 9 Be isotope ratios were measured by the 500 kV TANDY accelerator mass spectrometer at ETH Zürich and calibrated using the S2007 N 10 Be standard (Christl et al., 2013) (Table 3). The five chert samples from the Sfakia fan were cleaned by etching with hydrochloric and nitric acid and measured at the Cologne AMS facility (Dewald et al., 2013) relative to standards KN01-6-2 and KN01-5-3. 3.5 | Post-burial production and (paleo-) denudation rate calculations We calculated (paleo-)denudation rates with the CRONUScalc v2.1 tools (Marrero et al., 2016a). However, we modified the scaling factor functions (scalefacs1026, scalefacs36) to calculate scaling factors based on a pixel-by-pixel approach for the entire catchment area.
Catchment scaling factors were calculated on a 30-m SRTM DEM for modern and paleo-denudation rates. Time-dependent production rates were calculated with the Stone (2000) scaling and accounting for geomagnetic variations after Nishiizumi et al. (1989). We assume a bulk density of 2.65 and 2.2 g/cm 3 for the bedrock and fan deposits (Rodés et al., 2011), respectively, with an uncertainty of 0.1 g/cm 3 .
We developed a Monte Carlo routine that utilizes the CRONUScalc v2.1 production rate functions to estimate post-burial production ('PostPro') (Ott, 2022). First, the age of the sampled layer and age constraints on layers above are used to build a time-depth history, assuming that aggradation between two layers of known age The computation of post-burial production requires an age for the end of deposition above the sample. In most sampled locations, we used the mean age of the youngest deposited layer and the oldest incision constraint, which often corresponds to the depositional age of the next younger unit. However, at the Sfakia fan and a Stomio River terrace (Elafonisi), we also report an alternative scenario based on less conservative assumptions of surface abandonment. The analysis presented in the main text is based on the more conservative estimates for surface abandonment ages listed in Table 4. We report a detailed description of constraints on surface abandonment and uncertainty calculation in Table S2.  This marine formation consists of unsorted cobbles, pebbles, and sand derived from the PQ unit, mixed with algal constructions in a matrixdominated, well-cemented carbonate rock (Figure 3c). We interpret these deposits as the subaqueous part of the Stomio River delta system that has been uplifted above sea level and refer to this unit as the Stomio Formation. Chrissoskalitissa is constrained by a radiocarbon age from a snail shell within the channel that incised the fan at 11.3-11.2 ka (2σ range) ( Figure 3g). Two paleoshorelines near the Stomio River mouth, with

New chronological constraints
ISAs of 53 and 23 m, were dated to 58 AE 6 and 51 AE 8 ka, respectively.

Synthesis of new and published chronological data
Our new geochronologic ages and published observations constrain the timing of aggradation and incision at Elafonisi (Figure 4). Our luminescence age of Qt2 suggests deposition around MIS 5a/b. A paleoshoreline forming a paleocliff in Qt2 south of the Stomio River mouth was dated to 71 AE 7 ka (Ott et al., 2019b) and potentially marked the incision into Qt2. However, we dated a paleoshoreline at the same elevation to 58 AE 6 ka north of the Stomio River mouth. Both ages overlap within uncertainty, yet paleoshorelines on Crete have been shown to form during sea-level highstands and record relatively steady uplift rates (Gallen et al., 2014;Ott et al., 2019b;Robertson et al., 2019). Therefore, we assume the 71 AE 7 ka age lies closer to the true age because it matches the MIS 5a highstand, the most commonly preserved highstand on Crete (Ott et al., 2019b), and infer an uplift rate almost identical to the one derived from the paleoshoreline below (see Section 4.1.2). We note that the 58 AE 6 ka sample was collected at a shallow depth below the surface soil and younger grains may have been worked into the sediment from bioturbation (see Figure S1, sample USU-485). A paleoshoreline dated to 51 AE 8 ka and inset against Qt3 serves as a minimum age constraint for the end of

| Sfakia fan depositional chronology
As mentioned in Section 2, published geochronologic ages from the Sfakia fan show Qt1 deposition during MIS 5 and, after a short period of incision, Qt2 deposition from MIS 4 to 2 (Pope et al., 2008(Pope et al., , 2016 ( Figure 6). Currently, both units are being incised by the active channel, with up to 45 m of incision into Qt2 (Figures 7a and c). Incision commenced between 5.5 and 11.2 ka (Pope et al., 2008), and we dated nine carbonate cement rinds from the uppermost distal section of the fan (Figures 6 and 7d) to improve these estimates. The calibrated age ranges show a large spread from 3.7 to 11.1 ka ( Table 2). The age spread is likely related to multiple periods of cement formation (Geyh & Eitel, 1997). Moreover, the 'hard-water effect' that describes the incorporation of 'dead' carbon from local bedrock could bias the ages by up to one half-life (Shotton, 1972).
Therefore, we use the oldest cement ages to estimate fan surface abandonment in our alternative burial scenario (supplement) but apply

| Paleo-denudation rates
We estimate post-burial production for all cosmogenic nuclide samples ( Figure 8, Table 3, Figures S2 and S3) to calculate paleodenudation rates. For most samples, our estimated post-burial production is <30% of the total concentration; only for sample Sf-919-12, we estimate that 44% of the measured 10 Be and nearly half of the measured 36 Cl concentration can be explained by production after deposition ( Figure 8, Table 4). For samples El-919-3 and El-919-4 from the Chrissoskalitissa fan, the post-burial 36 Cl production is probably high. However, the uncertainties are large due to the high concentration ($300 ppm) of natural Cl. The calculation of post-burial production with our alternative burial scenarios shows that the difference in the final concentrations is <10% (Table S3).
We calculated (paleo-)denudation rates from our post-burialcorrected concentrations (Table 4, Figure 9). Denudation rates range from 108 to 1421 mm/ka, excluding the Chrissoskalitissa sample that yielded a concentration too low for a denudation rate calculation. At the Stomio River in the Elafonisi fan system, the modern denudation rate is 108 AE 10 mm/ka. Paleo-denudation rates are higher (136 AE 27 and 159 AE 44 mm/ka) but overlap at one-sigma uncertainty. The paleo-denudation rate from MIS 4 at the Livadia fan system is higher (343 AE 116 mm/ka), but from a significantly smaller catchment (Figure 2). At the Chrissoskalitissa fan, high concentrations of natural Cl cause large uncertainties in the post-burial production. The post-burial-corrected concentration of sample El-919-5 was too low for a denudation rate calculation. 36 Cl-derived paleo-denudation rates appear to be higher during MIS 2 compared to modern rates, yet due to the larger uncertainty in production rates, the modern and 23 ka paleo-denudation rates overlap within uncertainty.
At the Sfakia fan, most paleo-denudation rates are similar to modern rates, with some notable exceptions. The denudation rates from 10 Be and 36 Cl measurements agree for the modern stream sediment (264 AE 24 and 236 AE 20 mm/ka, respectively). However, three out of four 10 Be paleo-denudation rates are offset to lower values than the respective 36 Cl rates (Figure 9). 10 Be-derived denudation rates are similar to or slower than modern rates during MIS 2 and were faster at the MIS 3/4 boundary. 36 Cl-derived denudation rates show the same pattern but are offset to larger values.
Nuclide accumulation during denudation predicts a 36 Cl/ 10 Be ratio that is dependent on the denudation rate and the sample-specific chemistry. We calculated the predicted ratios for a range of plausible denudation rates (50-500 mm/ka) and production parameter uncertainties ( Figure 10). Apart from the modern sample, all samples exhibit 36 Cl/ 10 Be ratios below the expected values. A prolonged deep sample burial can explain ratios larger than expected due to the greater percentage of muon production for 36 Cl. Alternatively, long and shallow (20-200 cm) sample material burial can reduce the ratio ($3.5 to 6 depending on sample composition). However, our independent constraints on aggradation and incision do not allow for significant burial at depths shallow enough to significantly reduce the 36 Cl/ 10 Be ratios.
Another explanation for a low 36 Cl/ 10 Be ratio is the incomplete removal of meteoric 10 Be from our chert samples. It has previously been shown that adsorbed 10 Be can be challenging to remove from chert and lead to an excess of 10 Be (Zerathe et al., 2013). Without any other reasonable explanation, we assume that the incomplete removal of adsorbed 10 Be from the amorphous silicates led to the low 36 Cl/ 10 Be ratios. Therefore, we focus on the 36 Cl-derived denudation rates in our interpretation of the Sfakia fan evolution.

| DISCUSSION
We investigated the aggradation and incision timing and past varia- T A B L E 4 Post-burial production estimates and paleo-denudation rates. The numbers in the first column correspond to the location numbers in Figures 2 and 5 No. Sample Model surface age (ka) Post-burial estimate (at/g) Corrected 10 Be (at/g) Erosion rate (mm/ka) 10 Be F I G U R E 9 Post-burial corrected concentrations and (paleo-)denudation rates from Elafonisi and Sfakia fans. L = Livadia fan. One of the denudation rates in panel (b) from the Chrissisoskalitissa fan could not be calculated because the concentration was too low. Therefore, the denudation rate is assumed to be high, as indicated by the arrow. Samples from modern stream sediment are indicated with an asterisk (*). [Color figure can be viewed at wileyonlinelibrary.com] the last glacial-interglacial climate cycle. Below, we first discuss the controls on aggradation and incision in our alluvial systems before discussing our paleo-denudation rates. Subsequently, we examine regional changes in climate, sea level, and vegetation throughout the same period and discuss the dominant factors driving changes within these alluvial systems.

| Timing of aggradation and incision
Our geochronologic data show that aggradation at the Elafonisi fan sequence occurred during MIS 5b/a (Qt2), MIS 4 to early MIS 3 (Qt3), and MIS 2 (Qt4) (Figure 4). At the Stomio River mouth, no Qt4 deposit is preserved. Due to its large catchment area, we assume that the Stomio River was able to incise substantially during the low sea level of MIS 2 and likely deposited an MIS 2 fan that now lies offshore. This scenario is supported by the embayment of the coastline directly at the Stomio River mouth (Figure 2), which suggests river incision during eustatic regression. We did not date Qt1; however, the steep fansurface gradient suggests deposition during low sea levels. Based on its gradient and the other aggradation events broadly coinciding with the end of stadials, we assume the deposition of Qt1 during MIS 6. Incision events at the Elafonisi fan broadly coincide with interstadials and interglacial of MIS 1, 3, and 5a.
In contrast to the Elafonisi fan sequence, the Sfakia fan aggraded during most of the past glacial cycle (Pope et al., 2016). However, changes in the aggradation rate may still exist that we can derive from the dated sedimentary cross-sections published by Pope et al. (2008Pope et al. ( , 2016 (Figure 11). Combined with accumulation rates from the Chrissoskalitissa fan, dated in this study, aggradation was faster during late MIS 5 and 4 (90-60 ka), and especially during MIS 2 (29-14 ka) ( Figure 11a). During MIS 1, all fan sections are being incised. However, we note that apparent aggradation rates are typically affected by the timespan of measurement (Sadler, 1981). The timespans represented by our aggradation rates span one order of magnitude and show a moderate decrease with increasing timescale. We thus corrected the aggradation rates for a potential 'Sadler effect' by fitting a power law through the time scale-aggradation rate data (e.g. Gallen et al., 2015) (Figure 11b). Subsequently, we used this power-law fit to correct the aggradation rates in Figure 11a. For every aggradation measurement, we calculated an expected aggradation rate based on the power-law fit and subtracted the expected aggradation rate from the measured one (Figure 11c), similar to the approach of Kemp et al. (2020). These timescale-corrected relative aggradation rates show fast aggradation during MIS 4 and 2 and low rates during early MIS 5 and 3.
We thus conclude that at both study sites, aggradation increased during or towards the end of stadials (MIS 2,4,6), and incision or a reduction in aggradation rate occurred during interstadials and interglacials (MIS 1, 3, 5). Additionally, we compiled all geochronologic ages from Pleistocene alluvial fans and river terraces on Crete (n = 49, Figure 12a). The summed probability density distribution indicates two distinct peaks of aggradation ages, which correspond to MIS 2 and 4 ( Figure 12a) and align well with the timing of increased aggradation at the study sites.

| Paleo-denudation rates
Our paleo-denudation rates are slightly elevated compared to modern denudation rates, although most differences are within the margin of error ( Figure 9). The elevated uncertainty of the reported paleodenudation rates is partly due to the uncertainty in post-burial nuclide accumulation. However, the largest uncertainties arise from the high [Cl] content of the three samples at Chrissoskalitissa, where large uncertainty on thermal and epithermal neutron production rates complicates the interpretation of changes in denudation rate. In most other samples, the post-burial production component is minor, and even the two depositional scenarios we examined for several samples show a <10% difference in paleo-denudation rates. Therefore, we assume our paleo-denudation rate estimates to be robust.
At the Elafonisi fan system, paleo-denudation rates along the Stomio River are within the error of the modern rates measured. The denudation rates in the small Livadia and Chrissoskalitissa fan catchments ($1.5 km 2 drainage area each) are greater by a factor of 3-9 compared to the Stomio River (35 km 2 ). The Chrissoskalitissa catchment is steeper (32 ) than the Livadia (25 ) and Stomio (26 ) catchments and could explain why this location yields the fastest reported denudation rate so far for Crete (Ott et al., 2019a). However, both small catchments yield fast denudation rates, suggesting that the catchment area plays a role in the measured concentrations. In particular, the low paleo-denudation rates at Chrissoskalitissa could be explained by stochastic, low-concentration sediment input from mass wasting (Niemi et al., 2005). The Chrissoskalitissa catchment contains significant gypsum bedrock, which could promote slope failures.
Aggradation due to landsliding has been shown to occur in a similar carbonate catchment on Crete (Bruni et al., 2021). However, CRN denudation rates measured elsewhere have been demonstrated to be accurate and consistent, despite significant landsliding within

| Linking aggradation, incision, and denudation to base level, climate, and vegetation
Periods of pronounced aggradation align with stadials-sea-level lowstands-during the last glacial cycle (MIS 2, 4, 5b). Low sea level during these periods should favour incision in response to a base-level fall and seaward shift of alluvial depocentres. The brief period of incision observed at both study sites at the MIS 5a/4 transition might be an example of such base-level controlled entrenchment. However, F I G U R E 1 1 Accumulation rates of stratigraphic sections published in Pope et al. (2008Pope et al. ( , 2016 (Langgut et al., 2011;Wulf et al., 2018). We hypothesize that such variations in climate are the main control on the aggradation and incision behaviour of the studied fans and dominate the influence of base-level variations.
Aggradation during stadials and incision during interstadials has been observed in other Mediterranean regions (Macklin et al., 2002;Pope & Wilkinson, 2005;Wegmann & Pazzaglia, 2009;Zembo et al., 2009). For instance, Wegmann and Pazzaglia (2009) found that depending on the hillslope-derived sediment supply in the Northern Apennines, river aggradation or strath terrace formation occurred F I G U R E 1 2 Summary of aggradation and paleo-denudation data. (a) Individual (blue) and summed probability density functions of compiled geochronologic ages from alluvial fans on Crete (this study; Gallen et al., 2014;Holcomb et al., 2021;Pope et al., 2008Pope et al., , 2016. (b) Mean relative sediment accumulation rate from Figure 11C. (c) Aggradation and incision constraints from Elafonisi. (d) Paleo-denudation rates from this study. The Chrissoskalitissa (orange) denudation rates are scaled to the right y-axis for visibility. Rates from the Sfakia fan (light blue) and Stomio River (dark blue) are scaled to the left axis. Thick lines behind the dots indicate the approximate integration time of the cosmogenic nuclide sample. Dashed lines correspond to modern reference values. (e) Percentage of total tree pollen from Tenaghi-Philippon (TP, blue) in the northern Aegean (Wulf et al., 2018) and Ioannina (green) in northwest Greece (Roucoux et al., 2011). We calculated our age model for the Ioannina data based on linear interpolation between dated layers. Hence, the Ioannina pollen record's age uncertainty is greater than Tenaghi-Philippon. (f) Sea-level curve by Lambeck and Chappell (2001). [Color figure can be viewed at wileyonlinelibrary.com] during stadials with subsequent incision during the transition to interstadials. A compilation of aggradation periods in the Mediterranean by Macklin et al. (2002) documents major aggradation during MIS 6, the MIS 5b/a transition, MIS 4, MIS 2, and universal incision during MIS 1 at all studied sites. These aggradation periods align with the aggradation at the Elafonsi fan system and periods of increased aggradation at the Sfakia fan. About half of the sites in the Macklin et al. (2002) compilation lie inland. Therefore, the temporal correlation of aggradation and incision in the Eastern Mediterranean is additional evidence that both are mostly governed by climate-related changes in sediment transport, with only second-order controls by tectonics and sea level.
The aggradation behaviour of rivers on orbital time scales is a function of sediment supply and runoff-controlled transport capacity.
Our paleo-denudation rates show that hillslope denudation remained fairly constant throughout the glacial cycle, suggesting no major changes in hillslope denudation during the variations in aggradation and incision. This observation is in agreement with findings from the Northern Apennines, Central Europe, and the Himalayas, where Pleistocene denudation rates remain relatively stable (Cyr & Granger, 2008;Kapannusch et al., 2020;Schaller et al., 2002), despite shifts in river aggradation and incision (Wegmann & Pazzaglia, 2009).
However, the fluctuations in precipitation observed by paleoclimate studies in Greece suggest major variations in runoff-controlled stream transport capacity during the same time. We thus hypothesize that increased transport capacity due to wet conditions during interstadials and interglacials is probably the main driver of stream incision, whereas reduced transport capacity during stadials led to aggradation.
This interpretation suggests that streams in western Crete are transport-limited and sensitive to hydroclimate fluctuations.
Climate changes also led to large fluctuations in vegetation, which may impact sediment supply. Pollen records from northern Greece ( Figure 12e) and Israel show forest dominance during interstadials and interglacials and a switch to steppe conditions during stadials (Langgut et al., 2011;Margari et al., 2009;Roucoux et al., 2011;Wulf et al., 2018). Vegetation changes could have modulated the sediment supply to the channel by stabilizing hillslope material (Amundson et al., 2015;Torres Acosta et al., 2015). Therefore, shifts from forest to steppe conditions during stadials may have released pulses of hillslope sediment, leading to aggradation, similar to what can be observed in numerical simulations by Tucker and Slingerland (1997).
Faster hillslope transport due to less vegetation would decrease the measured cosmogenic nuclide concentrations (Anderson, 2015). However, this effect may be buffered by (1) the initial release of higherconcentration material that was stabilized, (2) sediment mixing within the catchment, which may dilute the input of low-concentration material from the hillslopes, and (3) the sluggish decrease in CRN concentration due to the integration time of CRNs. Our more variable paleo-denudation rates from the small Chrissoskalitisssa fan could indicate that small catchments with high denudation rates may be capable of registering a short-lived increase in denudation. In the larger Stomio River and Sfakia catchments, sediment storage is limited and integration times of cosmogenic nuclides are short (Figure 12d).
Yet, the addition of various buffering effects may be sufficient to damp a signal in cosmogenic nuclide-derived denudation rates to within uncertainty of modern rates, as observed in this study. Hence, it is likely that the large shifts in vegetation contributed to aggradation and incision by releasing and withholding hillslope sediment, respectively; however, the uncertainties and limited number of our paleo-denudation rates do not allow us to infer if changes in hillslope denudation played a secondary role in the aggradation behaviour of streams and alluvial fans.
Relatively steady hillslope denudation during aggradation-incision events has also been observed in the Himalayas (Kapannusch et al., 2020). A numerical model of 1-D alluvial river evolution by Simpson and Castelltort (2012) showed that variations in water discharge get amplified and transmitted downstream, leading to periods of pronounced aggradation and incision. In contrast, fluctuations in sediment discharge get damped on their way from source to sink. Our observations agree well with these findings. Observed variations in sediment input on Crete were minor and, if present, appear to have been damped over a length scale of a few kilometres. However, variations in hydroclimate left a strong imprint on the studied coastal alluvial systems, even outcompeting high-amplitude variations in sea level. The consistency between our data, theoretical expectations, and a few other studies suggests that alluvial sequences may rather record variations in hydroclimate than sediment input in many environments.

| SUMMARY AND CONCLUSIONS
Quaternary alluvial deposits at the Elafonisi fan complex, Crete, Greece, were deposited during MIS 2, 4, 5a/b, and likely also MIS 6. Similarly, aggradation rates at the Sfakia fan suggest the fastest deposition during MIS 2 and 4, with a brief period of minor incision at the MIS 5a/4 boundary and major incision during MIS 1. Despite large variations in aggradation and incision, millennial-scale denudation rates remained relatively constant in western Crete during the Late Pleistocene, except for one small catchment with elevated rates during MIS 2.
Periods of aggradation coincide with dry climate and steppe vegetation. We therefore attribute temporal changes in aggradation primarily to decreased runoff, potentially supported by a transient increase in hillslope sediment release due to a decline in vegetation cover. Major incision during MIS 1 is likely a response to increased rainfall and stream-transport capacity alongside denser vegetation stabilizing hillslope material. An early period of incision at the MIS 5a/4 transition likely occurred in response to rapid sea-level fall, suggesting that base level-driven incision events are superposed onto the climatically driven sequence of aggradation and incision.
Hence, we find that streams in western Crete show transportlimited behaviour, where variations in hydroclimate primarily modulate aggradation and incision while millennial-scale denudation rates remain relatively constant. In accordance with numerical models, we speculate that sediment-supply variations are only recorded in small and rapidly eroding basins, while changes in runoff dominate the alluvial record.