of Stochastic alluvial fan and terrace formation triggered by a high-magnitude Holocene landslide in the Klados Gorge, Crete

Abstract. Alluvial fan and terrace formation is traditionally interpreted as a
fluvial system response to Quaternary climate oscillations under the
backdrop of slow and steady tectonic activity. However, several recent
studies challenge this conventional wisdom, showing that such landforms can
evolve rapidly as a geomorphic system responds to catastrophic and
stochastic events, like large-magnitude mass wasting. Here, we contribute to
this topic through a detailed field, geochronological, and numerical
modelling investigation of thick (>50 m) alluvial sequences in
the Klados catchment in southwestern Crete, Greece. The Klados River
catchment lies in a Mediterranean climate, is largely floored by carbonate
bedrock, and is characterised by well-preserved alluvial terraces and inset
fans at the river mouth that exceed the volumes of alluvial deposits in
neighbouring catchments of similar size. Previous studies interpreted the
genesis and evolution of these deposits to result from a combination of
Pleistocene sea-level variation and the region's long-term tectonic
activity. We show that the >20 m thick lower fan unit,
previously thought to be late Pleistocene in age, unconformably buries a
paleoshoreline uplifted in the first centuries CE, placing the depositional
age of this unit firmly in the late Holocene. The depositional timing is
supported by seven new radiocarbon dates that indicate middle to late Holocene
ages for the entire fan and terrace sequence. Furthermore, we report new
evidence of a previously unidentified valley-filling landslide deposit that
is locally 100 m above the modern stream elevation, and based on
cross-cutting relationships, it predates the alluvial sequence. Observations
indicate the highly erodible landslide deposit as the source of the alluvial
fill sediment. We identify the likely landslide detachment area as a large
rockfall scar at the steepened head of the catchment. A landslide volume of
9.08×107 m3 is estimated based on volume reconstructions of the
mapped landslide deposit and the inferred scar location. We utilise
landslide runout modelling to validate the hypothesis that a high-magnitude
rockfall would pulverise and send material downstream, filling the valley up
to ∼100 m. This partial liquefaction is required for the
rockfall to form a landslide body of the extent observed in the valley and
is consistent with the sedimentological characteristics of the landslide
deposit. Based on the new age control and the identification of the
landslide deposit, we hypothesise that the rapid post-landslide aggradation
and incision cycles of the alluvial deposits are not linked to long-term
tectonic uplift or climate variations but rather stochastic events such as
mobilisation of sediment in large earthquakes, storm events, or ephemeral
blockage in the valley's narrow reaches. The Klados case study represents a
model environment for how stochastically driven events can mimic
climate-induced sedimentary archives and lead to deposition of thick alluvial
sequences within hundreds to thousands of years, and it illustrates the
ultrasensitivity of mountainous catchments to external perturbations after
catastrophic events.



12
friction angle. The Voellmy rheology is defined by Eq. (2), where μ is the frictional coefficient (equivalent to 13 ), ρ is the material density in kg m -3 , g the gravitational acceleration in m s -2 , v is the depth-averaged 14 flow velocity in m s -1 , and ξ is the turbulence coefficient in m s -2 . In short, adding to the basic frictional rheology 15 equation, Voellmy rheology includes a "turbulent term" which is dependent on flow velocity and the density of 16 the material and summarises the velocity-dependent factors of flow resistance (Hungr and Evans, 1996).

19
The input parameters for both rheologies need to be defined through back-analysis. Constraints on the parameters 20 may be deduced from the deposit's extent, the runout topography, and the material exposed along the sliding path.

21
Additionally, previous studies provide first estimates of reasonable input parameters in similar environments. A 22 common issue when modelling with DAN3D is that fluid pressure induces lateral spreading of a flow-like rock 23 mass already in the source area (Aaron and Hungr, 2016b). However, it is more reasonable to assume the rock 24 mass slides without much internal deformation in the rock avalanche's early stages. Therefore, a modified dynamic     (Table 3). Nevertheless, within these model runs, we managed to implement several parameter 53 combinations, and at least one resulted in a realistic landslide runout. Thus, the feasibility of the proposed 54 processes was tested successfully, and we achieved the primary goal of runout modelling.    Sharp geomorphic features, such as earthquake fault scarps and river terrace risers, smooth over time due to active 112 hillslope processes such as rain splash, creep, and tree throw (i.e., Carson and Kirkby, 1972;Selby and Hodder, 113 1993). This process can be approximated by simple diffusion models, which, when calibrated, have been used to 114 infer the timing of an earthquake or abandonment of a terrace (i.e., Culling, 1960;Fernandes and Dietrich, 1997).

115
Here, we use a linear diffusion model, the current and reconstructed initial morphology of the paleo-sea cliff

121
We construct a structure from a motion (SfM) digital surface model (DSM) for the Klados coastal fan sequence 122 from images collected from a DJI Phantom 3 drone using the AgiSoft photogrammetry software. The SfM point 123 cloud was converted to a DSM with a horizontal resolution of 15 cm (Fig. S5a, b). From the DSM, we extracted 124 two 10 m wide swath profiles from the modern shoreline and onto the upper fan surface in transection with 125 minimal vegetation (Fig. S5a, b). To minimize the impact of vegetation, we generate profiles using the minimum 126 values in the swath profile (grey lines in Fig. S5c, d). These profiles were used in the diffusion modelling.

127
We reconstructed the initial paleo-sea cliff morphology by performing a linear regression on the lower fan tread

142
The best-fit results for a given model run time all fit the observed profile equally well, and for this reason, we 143 only show one best-fitting profile (green line in Fig. S5c, d). There are some mismatches between the observed 144 and modelled profiles that are likely explained by processes not included in the simply linear diffusion model, but 145 despite these discrepancies, the results provide a reasonable approximation of the modern paleo-seacliff 146 morphology. The diffusion coefficients for our preferred Holocene (5 -7 kyr) abandonment timing of the seacliff 10 range between ~35 and 100 cm2/yr and those from Pleistocene (35 -45 kyr) abandonment timing range between 148 ~6 and 14 cm2/yr (Fig. S5e, f).

150
To assess if the best-fit results are reasonable, we compare them to a global compilation of diffusion coefficients