Articles | Volume 9, issue 4
Research article
17 Aug 2021
Research article |  | 17 Aug 2021

Controls on the grain size distribution of landslides in Taiwan: the influence of drop height, scar depth and bedrock strength

Odin Marc, Jens M. Turowski, and Patrick Meunier

The size of grains delivered to rivers by hillslope processes is thought to be a key factor controlling sediment transport, long-term erosion and the information recorded in sedimentary archives. Recently, models have been developed to estimate the grain size distribution produced in soil, but these models may not apply to active orogens where high erosion rates on hillslopes are driven by landsliding. To date, relatively few studies have focused on landslide grain size distributions. Here, we present grain size distributions (GSDs) obtained by grid-by-number sampling on 17 recent landslide deposits in Taiwan, and we compare these GSDs to the geometrical and physical properties of the landslides, such as their width, area, rock type, drop height and estimated scar depth. All slides occurred in slightly metamorphosed sedimentary units, except two, which occurred in younger unmetamorphosed shales, with a rock strength that is expected to be 3–10 times weaker than their metamorphosed counterparts. For 11 landslides, we did not observe substantial spatial variations in the GSD over the deposit. However, four landslides displayed a strong grain size segregation on their deposit, with the overall GSD of the downslope toe sectors being 3–10 times coarser than apex sectors. In three cases, we could also measure the GSD inside incised sectors of the landslides deposits, which presented percentiles that were 3–10 times finer than the surface of the deposit. Both observations could be due to either kinetic sieving or deposit reworking after the landslide failure, but we cannot explain why only some deposits had strong segregation. Averaging this spatial variability, we found the median grain size of the deposits to be strongly negatively correlated with drop height, scar width and depth. However, previous work suggests that regolith particles and bedrock blocks should coarsen with increasing depth, which is the inverse of our observations. Accounting for a model of regolith coarsening with depth, we found that the ratio of the estimated original bedrock block size to the deposit median grain size (D50) of the deposit was proportional to the potential energy of the landslide normalized to its bedrock strength. Thus, the studied landslides agree well with a published, simple fragmentation model, even if that model was calibrated on rock avalanches with larger volume and stronger bedrock than those featured in our dataset. Therefore, this scaling may serve for future modeling of grain size transfer from hillslopes to rivers, with the aim to better understanding landslide sediment evacuation and coupling to river erosional dynamics.

1 Introduction

Grain size is an essential parameter for understanding sediment transport and associated processes in river evolution or in hazards related to sediment pulses. For geomorphologists, it is increasingly considered to be an important parameter for the long-term incision of bedrock streams (Sklar and Dietrich2001; Cook et al.2013, 2014; Turowski2018), and it is an essential part of the sedimentological signal which is ultimately archived in stratigraphy (e.g., Armitage et al.2011).

Nevertheless, there are many processes that control the grain size distribution (GSD) delivered to rivers, and they are poorly understood (Allen et al.2015). In recent studies, models have been proposed that describe how weathering in the critical zone reduces the original size distribution of bedrock before the grains reach the surface (Marshall and Sklar2012; Riebe et al.2015; Sklar et al.2017). However, in active orogens with high erosion rates (>0.5 mm yr−1), landslides are likely the main providers of sediments to rivers (Hovius et al.1997; Struck et al.2015; Marc et al.2019), and a large fraction of sediment may reach the river only partially weathered. Indeed, the limits of models predicting soil GSD and the need to account for GSD derived from fractured bedrock was recently shown (Neely and DiBiase2020), although the role of mass wasting in delivering and further fragmenting bedrock particles was not explored. In those settings, understanding and modeling the controls on the landslide GSD should be an urgent goal, although, to date, it has only been addressed by a few studies. Indeed, in contrast to river sediments, for which many studies exist (e.g., Ibbeken1983; Whittaker et al.2011; Chung and Chang2013; Guerit et al.2014, 2018), landslide GSDs have rarely been measured; this is, in part, because the latter is considerably more difficult, time consuming and potentially dangerous than the former.

A few studies have measured and thoroughly discussed the GSDs of some large historical landslide or rock avalanches, often putting forward the various mechanisms of rock fragmentation and grain segregation to explain their data (see Crosta et al.2007, and references therein). Although interesting for their discussion in terms of rock mechanics, such case studies do not allow us to understand the regional variability in landslide GSDs nor to derive physical scalings that could pave the way to model the GSDs of material delivered to river networks by landslides. To our knowledge, only seven studies have reported detailed GSD measurements from multiple landslide deposits. A pioneering study reported the GSDs from 42 landslide dams across the Apennines, with a discussion on the methods used to derive the GSD but none on the controls of the GSD variability (Casagli et al.2003). Locat et al. (2006) presented GSDs from nine large (>100 Mm3) rock avalanches from Canada and the Alps, including various rock types, and analyzed these in terms of potential energy and fragmentation theories (Locat et al.2006). They found that the ratio of bedrock initial median block size, Di (estimated from fracture spacing), to the deposit median grain size, D50, was proportional to the change in potential energy per unit of volume, ρgH (where H  is the drop height of the center of mass, ρ is the rock density and g is the gravitational acceleration), normalized by the point load strength of the bedrock, σc, measured with a point load test performed on rock sample from the sites. Specifically their nine rock avalanches were best fit by a relation that could be recast as follows:

(1) D 50 = D i k 1 ρ g H σ c - k 2 ,

where k2=0.5 is an empirical threshold for fragmentation, and k1=83.3 is an empirical coefficient related to the conversion of potential energy into fragmentation energy and the effective breaking of particles. Thus, if the scaling is general for landslides, the deposit D50 should increase with Di and σc but decrease with H.

Subsequent studies often focused on the potential importance of landslide GSDs for understanding sediment transport dynamics and the expected GSDs at the outlet of basins, and they reported GSDs in Nepal, Japan, California and southern Italy (Attal and Lavé2006; Nishiguchi et al.2012; Attal et al.2015; Roda‐Boluda et al.2018). Several of them qualitatively underlined the factors influencing GSDs, such as the different lithological units (Attal and Lavé2006; Roda‐Boluda et al.2018) or the local hillslope gradient, as a control on the time spent in the weathering engine (Attal et al.2015). Recently, a study presented the GSDs of seven medium-sized rockfalls in Spain, showing that the bedrock block size and the deposit GSD could be related through a fractal fragmentation model (Ruiz-Carulla and Corominas2020). They found that potential energy was a main control on the fragmentation, but no clear correlation with rock strength measures emerged. They did not compare their model and results to the simple scaling proposed by Locat et al. (2006). Thus, none of these more recent works have attempted to frame the landslide GSD in terms of the competition between fragmentation energy and source rock strength, and the scaling for large rock avalanches has not been reproduced on smaller, more frequent landslides.

Based on these studies, we formulate two hypotheses: first, we suggest that Eq. (1) could be generalized to landslides of intermediate size and depth and, thus, that the landslide deposit D50 should increase with rock strength, σc, and the source material's median size, Di, but decrease with drop height, H; second, we hypothesize that materials mobilized by shallow landslides coarsen with the landslide scar thickness, T (i.e., Di increases with T), due to a reduction with depth of the fracture density of the bedrock (Clarke and Burbank2011) and/or of the degree of physical and chemical weathering experienced by particles (Cohen et al.2010; Anderson et al.2013; Sklar et al.2017). Testing these hypotheses seems essential to pave the way towards geomorphic models accounting for the GSDs of sediments transferred from hillslopes to rivers and from rivers to sedimentary basins (Allen et al.2015; Sklar et al.2017).

With these goals, and given the sparse amount of data on landslide GSDs, we performed detailed measurements on 17 recent landslide deposits in Taiwan. Taiwan is a prime example of an active mountain belt where landslides are the main supplier of sediment to rivers (Hovius et al.2000) and where reports of river GSDs exist in the literature (Chung and Chang2013; Lin et al.2014). Nevertheless, to our knowledge, comprehensive landslide GSD measurements are still lacking in Taiwan. Below, we report our measurements and discuss the source of variability in the GSDs within given landslides and across the whole dataset. We then discuss the validity of the two hypotheses stated above based on the GSDs of these landslides. We end by discussing the implications in terms of caveats and opportunities for GSD sampling and implications for fluvial sediment transport.

2 Data and methods

In this study, we report original GSDs for 17 landslide deposits from Taiwan (Fig. 1) as well as basic landslide information that we use to discuss controls on the GSD (Table 1). We also detail how we constrained landslide characteristics and measured GSDs for each deposit. Note that there are two landslides at the same site named LS-9o (for “old”) and LS-9n (for “new”), as the latter appears to have happened after the former (see Fig. 1d).

Figure 1(a) Hill-shaded elevation map of Taiwan, showing the main lithological units of the Central Range (based on the geological map of Taiwan available from the Taiwan Central Geological Survey) and the locations of the 17 sampled landslides deposits. (b–e) Pictures of some sampled landslides, where the yellow line is the approximate contour of the landslide (sometimes going beyond the pictures), and the dashed line indicates the transition from deposit to scar (it is only tentative when associated with “?”). In panels (b) and (e), the lead author is standing on the deposit for scale.

Table 1Landslide characteristics for the 17 surveyed deposits. Asterisks indicate landslides for which the geometry was estimated in the field rather than from satellite imagery. Due to its complex displacement, LS-12 has large uncertainties regarding the displacement of its center of mass. The abbreviations used in the table are as follows: Bsc – black schist (Tananao Formation), Sl/Sd – slate/sandstone (Lushan Formation), Sh/Sd – shale/sandstone (Nanchuang Formation), Msd – metasandstone (Shipachungshi Formation), Sh – shale (Chinshui/Cholan Formation), U – up, M – mid, T – toe, I – in, C – channel and n/a – not applicable.

Download Print Version | Download XLSX

2.1 Landslide characteristics

To quantify the variability in landslide GSD and its controls, we have targeted landslides with a known triggering date and landslides covering a broad range of areas (40 m2 to 0.1 km2) and lengths (10 to 400 m). Except for four small landslides (<1000 m2), which were opportunistically sampled close to larger neighboring ones, all landslides were targeted based on satellite imagery and chosen for the accessibility of their deposit.

Landslide type was difficult to assess, but most landslides could be called debris avalanches (Varnes1978), involving variable amounts of regolith and bedrock, although LS-13 and LS-14 could also be called rock falls. LS-12, the largest event, may rather be a deeper rock slump, with moderate displacement that is partly translational and partly rotational. Most landslides correspond to landslide polygons present in the Typhoon Morakot landslide inventory (Marc et al.2018); thus, they occurred in August 2009, about 5.5 years before they were surveyed in March 2015. Other more recent landslides were dated based on the time series of images available in Google Earth (see Table 1).

To assess variability in the GSD independent of rock type, 13 out of 17 landslides were chosen in the same geographic area, on both sides of the southern section of the Taiwanese Central Range, in relatively homogeneous lithological units composed of slate and slightly metamorphosed sandstone (Fig. 1, Table 1). LS-1 and LS-15 also occurred on moderately metamorphosed units, on both sides of the northern part of the Central Range, in black schist and in metasandstone intercalated with slate, respectively. The two remaining landslides both occurred in unmetamorphosed units made of alternating sandstone and shale for LS-10, which occurred in the emergent topography of Taiwan's southern tip, and in shales of the northwestern foothills for LS-16. In LS-16, many coarse rock fragments (>10 cm) were crumbling when touched, highlighting the weakness of this rock compared with the other units.

In an effort to constrain the mechanical strength of these units, we refer to measurements reported for 128 samples from the Chenyoulan catchment, both for the Nanchuang/Nankang Formation, which extends to the southern tip of Taiwan and contains LS-10, and for the metasedimentary units of the Shipachungshi Formation, where LS-15 occurred (Lin et al.2008). The unconfined compressive strength of Nanchuang sandstone ranged from 29 to 117 MPa (mean of 70 MPa), whereas the Shipachungshi metasandstone ranged from 45 to 179 MPa (mean of 100 MPa) (Lin et al.2008). However, in the Nanchuang Formation, sandstones alternate with weak shale (strength below 10 MPa), with an equal proportion of each. In contrast, the Shipachungshi metasandstone is intercalated with less frequent slates, often stronger than the Nanchuang shales (Lin et al.2008). These measurements clearly make the case for highly variable rock strength and are far from encompassing the potential diversity of the rock types sampled by the studied landslides. Our goal here is not primarily to constrain the rock strength of individual landslides but instead to estimate the relative strength of diverse units. Based on the measurements reported above, we make two assumptions: first, that the shales and sediments hosting LS-16 and LS-10 may be 7–13 times and 2–4 times weaker, respectively, than the metasediments hosting LS-15; and, second, that the slates and metasandstones in the Lushan Formation have similar strength to those in which LS-15, as well as LS-1, occurred and, thus, that these landslides can be compared without normalizing for strength.

Geometric landslide metrics were obtained from high-resolution satellite imagery available in Google Earth except for four deposits (LS-4, LS-9n, LS-13 and LS-14), which were too small to be clearly distinguished on the imagery, and which had their dimensions approximated from field observations only, using a laser rangefinder. Area was obtained by hand mapping the whole disturbed zone on the imagery. Length refers to the downslope length between the highest and lowest point of the polygon. The elevation difference between these two points, estimated from the elevation data of Google Earth (in Taiwan mostly the 30 m Digital Elevation Model, DEM, from the Shuttle Radar Topography Mission, SRTM, predating all of the studied landslides), defined the maximum drop height. Physically, the potential energy change in Eq. (1) is related to the displacement of the center of mass. Thus, we also estimate a length and drop height from the center of mass of the scar to that of the deposit, estimated from Google Earth imagery. An estimate of the pre-failure scar gradient could be derived from the scar's approximate length and height difference.

Beyond plan-view metrics, we must also estimate landslide scar volume in order to constrain the landslide scar depth. For a few landslides, the deposit volume could be approximated as a fraction (a quarter to a half) of a cone, for which a volume estimate could be obtained as πR2h/3, where R and h are the respective approximate radius and height of the cone that were estimated in the field. This simplified geometry was only suitable for LS-3, LS-4, LS-7, LS-9n, LS-10, LS-11, LS-13 and LS-15, and yielded only a first-order “field volume” estimate (Table 1). For the other deposits, we had to rely on scaling relationships between the scar area and volume, with the additional complexity that the lower extent of the scar area could not be clearly assessed in most cases, even with high-resolution imagery. Thus, assuming the mean scar aspect ratio from a global database applies to the surveyed landslides (Domej et al.2017), we estimated scar areas as As=1.5Ws2, where Ws is the scar width obtained by measuring the extent of the landslide in the direction orthogonal to flow, in the upper part of the failure only. With As, we estimated a maximum and minimum landslide volume using empirical scaling relationships of the form V=αAsγ, with different parameter values assuming the scar was mobilizing soil or bedrock, respectively. We used γ=1.262±0.009 and log10(α)=0.649±0.021 for landslides in soil, and γ=1.41±0.02 and log10(α)=0.63±0.06 for landslides in bedrock (Larsen et al.2010). We then derived the upper and lower estimates of landslide mean scar depth as the ratio of volume to scar area for bedrock and soil, respectively. Although approximate, this scaling is still preferable to using total landslide area, as it removes bias in the volume estimates associated with variable runout length, which is difficult to constrain from field or satellite observations (e.g., Marc et al.2019). Where available, the volumes estimated from the field mostly fall within the bracket of the volumes estimated from global scaling relationships (Fig. S1 in the Supplement), lending some support to this approach.

Nevertheless, for the deposits where we could not obtain a field estimate, a better constrained volume estimate could be obtained by choosing one of the two scaling relationships. We note that the field volume estimate of LS-3 and LS-4 is similar to the estimate from the soil scaling relationship (Fig. S1). This is consistent with the observation that they were composed of rock debris with a yellowish color that indicated advanced weathering, and contained fresh vegetation debris (see Fig. 1). For LS-7 and LS-11, which clearly involved mostly fresh bedrock, the field volume estimates match better with the bedrock scaling relationship (Fig. S1). Thus, where field volumes were lacking we used the bedrock estimate for the largest landslides (Ws>50) within which the rock looked mostly fresh (i.e., LS-1, LS-2, LS-8 and LS-12). Some other landslides (LS-5, LS-6, LS-9o and LS-16) featured a mixture of soil and rock material. Consequently, we used the average of the soil and bedrock scaling to estimate their volume. The best estimate for each landslide was divided by its scar area to obtain an estimate of scar thickness (Table 1).

2.2 Grain size count

GSDs were obtained using grid-by-number sampling, following established protocols developed for measuring riverine GSDs (see Kellerhals and Bray1971) and subsequently applied to landslide deposits (Casagli et al.2003; Attal and Lavé2006). We extended survey tapes along an elevation contour over a substantial portion of the deposit width (10 to 50 m) and sampled grains along the tape at a constant interval, recording the size bin of the b axis, measured with rulers. We used bins following a half-phi scale (power of 2 by 0.5 increments) with the smallest bin encompassing all grains finer than 2 mm. When grains could not be moved, we considered the smallest of the two visible axes as the b axis. The sampling interval was 0.5 m in most cases but was adjusted to 1 m for deposits where many meter-scale boulders where present (LS-2s, LS-13 and LS-14) to avoid having to count many grains several times. We then moved the line in parallel, upslope by 1 to a few meters depending on the deposit's dimensions and local topography, and repeated the count. Most slides were sampled with 6–10 survey lines, allowing us to cover a substantial fraction of the deposit (often 30 % to 60 %), with total counts typically including 200–400 individual grains. This approach also allowed us to sample different sections of the deposit when a spatial segregation was visible (LS-3, LS-8, LS-9n and LS-10) and to quantitatively assess this spatial variability in grain size (see Ruiz-Carulla et al.2015).

Figure 2(a) Range of D50 and D84 for various studies, excluding the coarsest and finest distributions of each study. (b, c) Cumulative distribution function for the 17 sampled landslide deposits. Note that we show the surface and inner distributions for LS-5, LS-8 and LS-2; thus, 20 distributions are shown in total. For visibility, some lines are dashed and the distributions are shown in two panels. Vertical lines are the approximate boundary for grain transport by suspension and bed load, for a flood associated with a fluid shear stress of 220 Pa (see Sect. 4.3).


Specifically, we observed and measured spatial variations in the GSD on the surface deposit for four landslides (LS-3, LS-8, LS-9n and LS-10). In addition to the surface GSD, we also measured the GSD of a section of the deposit likely to represent the interior of the deposit for three of the landslides (LS-2, LS-5 and LS-8). Below, we explain how we could independently measure various GSDs on these landslide deposits. For LS-8 and LS-2, we counted grains on the vertical banks of a 2 m deep erosional gully incising the deposit, and on a debris fan next to and below the road that had been cleared from the deposit, respectively. Thus, the former case allowed us to survey the internal GSD in place, whereas the latter likely represents a remixing from surface and internal parts and, thus, must be closer to the inner GSD than what would be derived from surface measurements only. Note that the only undisturbed deposit for LS-2 was the one in the transport channel, where a carapace (a layer of very coarse grains, Crosta et al. (2007)) seems to have formed (Fig. S2). Finally, on LS-5, we separately measured a debris fan and the terminal section of a channelized deposit which was visibly coarser (Fig. S3). In this case, given the age of the deposit and its direct contact with the floodplain, it is plausible that the deposit was partly eroded and that the fan may be a mixture of internal and superficial material, whereas the higher up channel section may be more representative of the original surface of the deposit. This will be further discussed in the following during the examination of segregation. In any case, to differentiate between the GSD considered to represent the interior or surface of the deposit, we add the letter “i” or “s” after the name, respectively (Fig. 2).

Figure 3Examples of three types of heterogeneity in grain size distributions within the same landslide deposit: (a) lithological difference within the deposit of LS-7 with sandstone grains coarser than slate grains; (b) downslope differences between the upper, middle and lower part of the deposits of LS-3 and LS-9; (c) difference between the surface and inner part of the deposits of LS-2, LS-5 and LS-8.


Additionally, on the deposit of LS-7, we could distinguish (by visual inspection) grains made of slate, which were dark, elongated and without visible internal structure, from grains made of metasandstone, which were lighter, more cubic and with visible internal grains. We counted them separately as we found them over the deposit. In many other deposits, a large majority of grains either looked fairly homogeneous or, due to time constraints, a systematic count within different rock types could not be done. To study the variability between various landslides, we obtained an overall surface GSD by summing the grain counts from both lithologies of LS-7 and from the different sectors of the deposits with spatial segregation. For these cases, we did not use an area-weighted sum (Ruiz-Carulla et al.2015) because the upper, middle (when differentiated) and lower sections of the deposits represented roughly similar proportions of the surface of the deposits, and we obtained count variations below 10 % from the different subsections (Table 1, Fig. 3). When we measured both internal and superficial GSDs, we had to select one of them as representative for the comparison to other landslide properties (see Sect. 3.2).

3 Results

3.1 Landslide grain size distributions and their internal variability

Before averaging the spatial variability, the landslide GSDs were found to have 50th and 84th percentiles ranging from about 15 to 200 mm and about 60 to 600 mm, respectively. This is consistent with the range of observations from previous studies, except the large rock avalanches from Locat et al. (2006) and the volcanic rock avalanches from Crosta et al. (2007), which were about 10 times coarser and finer, respectively, than all other studies. LS-2s and LS-16 are much coarser and finer than the rest of the studied landslides, respectively. Interquartile ratios vary between 3 and 15, but we note that 13 out of 20 GSDs have an interquartile ratio of 3 to 6, whereas only LS-1, LS-3, LS-5i, LS-5s, LS-8i, LS-15 and LS-16 have larger spreads (Fig. 2). All distributions seem unimodal except LS-16, which was found to have more than 40 % of the grains finer than 2 mm and likely contained a second, submillimetric mode that could not be constrained by our methods. Grain size distributions can often be well described by a Weibull or lognormal distribution (Ibbeken1983). For the studied landslides, eight GSDs are better fit (according to both Kolmogorov–Smirnov and Anderson–Darling statistics; Stephens1974) by a Weibull distribution (LS 2s, LS-3, LS-4, LS-5s, LS-5i, LS-6, LS-9n and LS-9o; Figs. 2b and S4), whereas all others are better fit by a lognormal distribution (Figs. 2a and S5). Note that LS-16 is poorly fit by both distributions. These two subgroups imply that we cannot prescribe a given distribution form to model landslide GSDs, and future work may aim at understanding why some landslide GSDs obey one or the other distribution. In any case, we refrain from using distribution parameters and will continue to discuss results based on empirical descriptors (i.e., median, interquartile ratio).

GSDs within a single landslide deposit were often heterogeneous: in one case, they were associated with differences between grains of different rock types (slate and sandstone in LS-7), whereas in seven other cases, they were associated with spatial variability across the landslide (Figs. 3 and S6). For LS-7, the slate pieces have grain sizes about three times smaller than the sandstone for a given quantile of the GSD, with a similar distribution shape. The slate grains were typically elongated platelets (i.e., with a axes about 3 times longer than the b axes, and c axes much smaller than the b axes), whereas the sandstone grains were cubic and slightly more abundant than the slate grains (N=196 vs. N=167). We observed (with the naked eye) downslope segregation, i.e., an increase in sediment coarseness from the apex to the toe of the deposits (Ruiz-Carulla et al.2015), in four cases. The strongest segregation was found to occur in deposits LS-3 and LS-9n, where the upper part of the deposits have grains 5-10 times finer than the lower part of the deposit, without substantially changing the shape of the distribution. Deposits LS-8 and LS-10 were found to exhibit a more subtle segregation, with the upper part of the deposits having distributions that were a factor of 1.5–2 finer than the toe of the deposits (Fig. S6). The toe of LS-10 also displayed D50 and D84 values that were twice as coarse as at its apex, consistent with other cases, but also had more fine grains, with about 10 % of grains finer than 2 mm vs. less than 5 % at the apex.

In two cases, we could separately measure the superficial and internal GSD. For LS-8, we observed that the superficial GSD had D16=20, D50=40 and D84=120 mm, whereas the internal GSD had D16=3, D50=10 and D84=50 mm. At the abovementioned location, the superficial deposits had almost no fine sediment below 2 mm, whereas the internal body had more than 10 % of fine sediments. Thus, the internal GSD had quantiles that were 10 to 20 times finer than the channel carapace – the largest difference observed in terms of internal variability. Note that the carapace also had a coarser GSD than any other measured landslide deposit in our study. In spite of this large difference, we note that the internal GSD still only had about 3 % of grains finer than 2 mm. These two examples clearly show that the superficial GSD can be substantially different from the internal GSD, both in terms of fine grains (<2 mm) and for coarse to very coarse grains (10 to 100 mm).

Last, in the case of LS-5, it is not entirely clear if the two distributions represent vertical segregation or superficial spatial variability. Given that the fan has a D16=4 and a D50=20 mm, about 3 times finer than in the channel, but an almost identical D84 of around 200 mm, we consider it to likely be an internal or mixed GSD.

3.2 Relationship with landslide properties

The percentiles of the GSDs are highly correlated with linear correlation coefficients of R2>0.9 between D50 and D16, D25, D75, D84 and D90 (Fig. S7). We note that more scatter is present for D16 and D90, suggesting that D50 is a good proxy for the bulk of the distributions but does not completely capture the variability in their tails. The interquartile ratio (here D75/D25), which characterizes the span of grain size in the distribution, is also independent of the other percentiles (Fig. S8).

For comparison with landslide properties, we used a single GSD for each landslide. For lithological variations or spatial variability on the surface, we averaged the various GSDs from each landslide. For the three landslides for which we have both an internal and superficial grain count, we identified a single GSD that was most relevant for comparison with the other deposits. For LS-8, we considered the GSD at the surface, to be consistent with all other cases. For LS-5, we considered the coarser distribution from the channel as more representative of the surface deposit. If the fan of LS-5 is representative of the deposit and if the channel is coarser because of some spatial segregation, the percentiles of its GSD would be about 2–3 times finer than the ones that we have selected (Fig. 3). In contrast, for LS-2, we considered the internal GSD, as the superficial measurement recorded only what seems to be a carapace that overrepresents coarse grains. Indeed, with about half of the distribution made of boulders (>0.5 m), we consider that segregation to be exceptional within our dataset. Thus, in using LS-2i to study inter-event variability, we assume that none of the other landslides had a carapace with similar strong sorting. Nevertheless, LS-2i percentiles may be biased towards finer grains, when compared with surface deposits of the other landslides, and we assume the bias could be up to a factor of 2–3, based on LS-5 and LS-8 (Fig. 3).

According to Eq. (1) and assuming similar initial regolith material for all landslides, D50 should decrease linearly with the ratio of drop height, H, to bedrock strength, σc. Indeed, log-transforming and fitting D50 against H for the 15 deposits from the metasedimentary units, which we expect to have relatively similar strength, we obtain R2=0.71 and a power law exponent of −0.64 (Fig. 4). Including LS-10 and LS-16 from weaker units yields a substantially poorer fit (R2=0.31). Thus, to account for their weaker strength, we rescaled these two landslide drop heights by a factor of 3 and 10, respectively, reflecting the central values estimated from strength measurements (see Sect. 2). We then obtained a correlation coefficient for all of the landslides (R2=0.71; N=17) and a best fit power law exponent of −0.78. For landslides in metasedimentary units, D50 is also negatively correlated (with a larger scatter R2=0.5–0.55) with landslide size metrics (area, width, volume and depth). However, we note that, for this dataset, these metrics are also strongly correlated with the drop height (R2=0.56–0.66; Fig. 5). Given that we would expect larger and deeper landslides to mobilize fresher and coarser grains and, thus, to have a positive correlation with D50, these negative correlations may simply reflect the fact that deeper landslides have a larger drop height in our dataset, and that the effect of drop height on decreasing D50 is more important than the effect of landslide size on sourcing coarser material.

Figure 4D50 for the 17 landslides (color coded by scar thickness) in this study against the drop height of their center of mass. Results are similar for maximal drop height. The best fit (solid black line) and R2 only consider the drop height rescaled by the landslide rock mass strength (black circles). The red circles show the original drop height for LS-10 and LS-16. Vertical bars show a factor of 2 uncertainty for LS-5s and LS-2i, for which there may be vertical segregation (see text). Horizontal error bars represent uncertainties in the drop height of LS-12 and the strength normalization of LS-10 and LS-16 (see Sect. 2).


Figure 5Correlation between the principal geometric dimensions of the surveyed landslides, drop height, scar width and scar depth. Note that scar depth is computed directly from scar width for all the landslides for which we do not have volume estimates from the field, i.e., all except LS 3, LS-4, LS-7, LS-9n, LS-10, LS-11, LS-13, LS-14 and LS-15 (see Sect. 2).


For the spread of the distribution, characterized by the interquartile ratio, we did not find any substantial correlation with any of the landslide variables. Even the rock type (or rock strength) does not seem to have an impact on the GSD spread, with several landslides in metasedimentary rocks with very large spreads (LS-3, LS-5 and LS-1), whereas the two landslides in unmetamorphosed units are on both ends of the spectrum. Thus, more data are needed to understand the spread of the landslide GSDs.

4 Discussion

In the following discussion, we propose that the variability in landslide D50 can be reconciled with the fragmentation scaling of Locat et al. (2006) (i.e., D50 decreases with the ratio of drop height to bedrock strength, as in Eq. 1), when accounting for regolith coarsening with depth (e.g., Cohen et al.2010). We then detail processes that may lead to grain size segregation within a given deposit as well as practical implications for sampling. Last, we explore the implication of the measured GSD on sediment transport and evacuation from river channels.

Figure 6(a) Grain size as a function of depth in the regolith inspired by the weathering model of Cohen et al. (2010), used to estimate the original median grain size, Di, mobilized by landslides with different thicknesses. (b) Reduction ratio, using the modeled Di from panel (a), against the potential energy normalized by point load strength estimate for the 17 landslide deposits in this study. The error bars are the same as in Fig. 4.


Before turning to these points, we recall that some of our variables that are important for modeling landslide D50 (Sect. 4.1) are only first-order estimates, if available at all. First, we rely on rough estimates of the landslide geometry (e.g., drop height, volume and depth). Better characterization in the field would have required more field work and was limited by accessibility or elaborate construction of DEM based on lidar or drone photogrammetry, which is difficult to perform and limited by the lack of accurate pre-failure DEMs. Thus, for the sake of this first study, we think that having consistent, first-order estimates of these metrics is sufficient to test the dependence of the D50 on landslide geometry. The difficulty of accessing many scars as well as the mixed origin (i.e., weathered regolith and bedrock) of several landslide sources also meant that we could not practically measure the source materials' median grain size, Di. Nevertheless, we propose below that variability in Di may be captured with existing weathering models.

4.1 The importance of fragmentation and source material initial grain size

Here, we discuss the hypothesis that Eq. (1), proposed and validated by Locat et al. (2006) for large rock avalanches, can also be used for smaller, shallower landslides made of a mixture of regolith and bedrock. For the 17 Taiwanese landslide in our study, we found that, within a given lithology, drop height seems to be a first-order control on the landslide deposit median grain size (Fig. 4). We also found that by rescaling the drop heights by their weaker rock strength, LS-10 and LS-16 were consistent with the trend defined by the stronger metamorphosed units. These observations qualitatively agree with Eq. (1); however, quantitatively, the best fit between H and D50 was not linear but a sublinear power law. Given that we observe that H and the landslide scar thickness, T, are correlated in our surveyed landslide (Fig. 5), this discrepancy with Eq. (1) could be resolved if Di, which we could not measure, increases with T. Models describing the size of particles in a soil or regolith predict upwards fining of grains from the bedrock to the surface due to an increase in the degree of both physical and chemical weathering (Cohen et al.2010; Anderson et al.2013; Sklar et al.2017). In bedrock, fracture density estimated from seismic wave refraction was also found to decrease nonlinearly from the surface to a depth of 5–10 m (Clarke and Burbank2011). Given that soils are often thin in Taiwan and represented a small proportion of the mobilized material, we consider that physical weathering is likely dominant. Here, we consider that the physical weathering rate (i.e., the rate of particle breakdown) can be modeled with an exponential decay from the surface, with a characteristic length scale of λ=2 m, consistent with previous modeling (Cohen et al.2010; Anderson et al.2013). Assuming the regolith median grain size at a depth zDr(z), to be proportional to the integral of the weathering rate, we modeled it as Dr=Db(1-exp(-z/λ)), where Db is the unweathered bedrock block size, producing a rapid variation near the surface consistent with published models for physical weathering (Cohen et al.2010; Anderson et al.2013). Averaging Dr from 0 to T, the mean scar thickness, and choosing a value for Db, we can obtain Di=Db(1-λ/T(1-e-T/λ)) for each landslide (Fig. 6a).

For comparison with the prediction of Locat et al. (2006) (i.e., Eq. 1), we have to assume a value for Db and use point load strength values (σc) of 0.5 MPa for LS-16, 1.5 MPa for LS-10 and 5 MPa for the other landslides in metasediments. These point load measurements are consistent with the typical unconfined compressive strength of intact rock (Lin et al.2008) after dividing them by 20 (Chau and Wong1996). With these strength values and assuming Db=1900 mm, the ratio of the modeled Di to the measured D50 agrees with Eq (1) (R2=0.9), even though its coefficients (k1 and k2) were calibrated on rock avalanches with rock strength and median grain size that were orders of magnitude larger than the ones in this study (Table 1, Figs. 2a and 6). All of the surveyed landslides are within a factor of 2 of the predictions from Eq. (1) – even LS-12, which likely had a different deformation style than the other landslides. This good agreement is not very sensitive to how we estimate landslide volume and scar thickness (Fig. S9). Note that by assuming Db=1900 mm we have matched the D50 of surface deposits, which may slightly overestimate the representative grain size relative to the whole landslide deposit, due to kinetic sieving or fine removal by surface runoff (Fig. 3c and discussion below). Based on LS-5 and LS-8, the inner D50 may be 2–3 times finer than its surface counterpart; thus, field measurement of regolith and bedrock GSD may need to be compared to a model with Db∼600–1000 mm, to which uncertainty in rock strength (i.e., probably a factor of 2) should also be added. Nevertheless, although uncertain, these latter values are intermediate between fracture measurement on surface outcrops in the US (Neely et al.2019; Verdian et al.2020), ranging from 10 to 400 mm, and the bedrock block size measured by Locat et al. (2006) on the scar of a large rock avalanche, ranging from 600 to more than 10 000 mm. Db∼600–1000 mm also matches quite well with the range of D50 to D75 (450 to 900 mm) found in the carapace of LS-2 and may be a fair first-order estimate of the original bedrock block size (Crosta et al.2007). We also note that LS-10 and LS-16, which occurred in weaker bedrock, may be expected to have a finer Db than the other slides. We are not able to constrain this, but we note that even with a Db 3 times finer, these two slides would only be a factor of 2 below the predicted D50 (Fig. 6).

We conclude this section by underlining that more measurements, especially of source rock block size and strength, are needed to fully demonstrate the applicability of the fragmentation theory presented by Locat et al. (2006). Nevertheless, we suggest that such fragmentation theory is applicable to understand and predict landslide GSD in a wide range of contexts, at least for rock, soil, and mixed avalanches and generally disrupted slides, which are the most commonly triggered (Keefer1984). Further, our observations suggest that Eq. (1) can be generalized to account for an exponential reduction in regolith grain size towards the surface (Cohen et al.2010; Anderson et al.2013), yielding the following:

(2) D 50 = 1 - λ T ( 1 - e - T / λ ) D b k 1 ρ g H σ c - k 2 ,

where Di has been replaced by a term depending on the “fresh” bedrock median size (Db), the length scale of weathering decay (λ) and the landslide thickness (T). In a sense, Eq. (2) supports previous qualitative statements on the importance of rock type (Attal and Lavé2006; Roda‐Boluda et al.2018), which may physically relate to rock strength and regolith block size. Additionally, Eq. (2) combines the concept of physical weathering with the process of fragmentation, controlled by drop height, the latter of which is less often considered in the geomorphological community. Future studies should also clarify whether Db is controlled by σc, which would make the equation more nonlinear but would reduce the number of parameters to constrain. More complex models of fragmentation have been used to predict landslide GSD (De Blasio and Crosta2014; Ruiz-Carulla and Corominas2020) and may be better suited to model the full GSD, but Eq. (2), provided it is further validated, opens various interesting perspectives. For example, it suggests that seismically triggered landslides, which occur more often near ridges than rainfall-triggered landslides (see Meunier et al.2008; Rault et al.2019) and are, thus, expected to have higher H, are more likely to deliver finer grains to the river systems, assuming they have a similar size and depth distributions (e.g., Marc et al.2019). More generally it highlights the need for an investigation on how geomorphic factors (e.g., hillslope height, steepness and shape) modulate landslide runout, drop height and connectivity to channels. Comparing hillslopes in various landscapes could be easily attempted based on comprehensive landslide inventories (Tanyaš et al.2017; Marc et al.2018). Equation (2) would also be well-suited for landscape scale modeling of the input of various grain sizes into rivers (e.g., Benda and Dunne1997; Carretier et al.2016; Neely and DiBiase2020) and, thus, for providing better coupling of landslide and river dynamics in landscape evolution models (Campforts et al.2020; Egholm et al.2013).

4.2 Controls on the internal variability in the GSD and implications for future sampling

We found three sources of internal variability in landslide GSD: one associated with the lithology of the individual grains, as reported for Himalayan landslides by Attal and Lavé (2006), and two related to the location of the grains on or in the deposit, as reported for various rock avalanches (Crosta et al.2007; Ruiz-Carulla et al.2015). We discuss these observations first in terms of implications for bias and sampling procedure, and second in terms of physical process causing them.

The lithological difference is not likely to be a bias as long as the grains of different lithologies are randomly distributed in the deposit: their sampling frequency should represent their relative abundance in the deposit. Spatial segregation on the surface of the deposit implies that – in order to ensure a representative GSD – the sampling method should ideally be performed across most of the deposit, or at least over the different subunits of the deposit, before doing a weighted average with their relative area of contribution. Measurement based on sieving at a single site or local grain counts along a line or over a fraction of the landslide area may misrepresent the GSD and should be avoided. For large deposits where access is difficult, the use of pictures from a drone may help with checking for segregation and potentially allow one to reproduce the grid-by-number counting method using image analysis (e.g., Casagli et al.2003; Attal and Lavé2006). However, this requires the scaling of each drone picture and, thus, the deployment of reference objects across the deposit, which is not always practical; moreover, such sampling will be unable to resolve fine grains (<30–100 mm). In contrast, in the presence of a vertical segregation, where superficial and inner GSDs differ, it may be very difficult to estimate a GSD that is representative for the whole deposit. When possible, targeting the banks of incised gullies may offer a good opportunity to characterize the subsurface of the deposit. Some applications mainly require the subsurface GSD – for example, modeling the weathering of freshly fragmented bedrock in the landslide deposit and how they can contribute to solute fluxes (Emberson et al.2016a, b). In contrast, the surface grains matter for sediment transport, and armoring may limit the mobilization of deeper finer grains. Additionally, in the case of a carapace, the question of how to combine the two end-member distributions would require an estimate of the relative thickness of the two end-member GSDs, which may be challenging. In the case of a less extreme segregation, as observed for LS-8 and probably LS-5, the proportion of coarse grains (>200 mm) was similar on the surface and inside the deposit, and only the medium and especially fine grains were more abundant inside the deposit.

The process of kinetic sieving (Savage and Lun1988; Gray2018) is expected to cause vertical segregation (i.e., a coarser surface and finer subsurface) in granular flows and a downslope segregation when shear is present, leading to boulder fronts as for LS-3. However, it should be noted that segregation is favored by transport along moderate slope gradients and tends to disappear for very steep chutes (Vallance and Savage2000). Although our gradient estimates are very rough, segregation mostly occurred for landslides with large transport distance, estimated as L2+H2, and the least steep slopes (Table 1, Fig. S10). This excludes LS-12, which likely has complex displacement, and LS-16, for which the weak and clay-rich lithology, prone to form agglomerates, may not behave like a typical granular material. Nevertheless, it seems hard to use kinetic sieving to explain why LS-9n was so clearly segregated downslope, in spite of its very modest size and displacement. Instead, we could hypothesize that the episodic reactivation of the scar and channel chute may have sprayed the deposit with finer debris on some landslide deposits, depositing preferentially near the apex of the deposit. Such a mechanism might have happened on most of the landslides that we have sampled (given their ages), but we cannot currently constrain its relevance without repeated monitoring of the deposits, which is left for future studies. Alternatively, for old deposits, it is likely that fine materials could have been washed away by repeated storm events. This progressive washing of the fine grains would be consistent with the fact that the superficial deposits are very poor in fine materials but have a proportion of coarse blocks fairly similar to the internal part of the deposit (e.g., LS-5 and LS-8; Fig. 3). In these two cases, kinetic sieving may have been limited (although likely present in LS-8 to explain some downslope coarsening) and fines may have been preferentially washed out. On various parts of the LS-11 deposit we did find finer materials when scraping off the top layer of gravels, consistent with this hypothesis. If such a process is expected to happen on all landslide deposits, superficial measurements of very fresh landslides may represent the bulk of the material (e.g., LS-15 and LS-3, as they most recently failed and have a high proportion of fine grains), and older deposits may require some correction, as medium to fine grains may be underrepresented.

To conclude this discussion, it seems clear that several physical processes can add complexity to landslide deposit GSDs and that deconvolving them and applying a process-based correction is not straightforward. More datasets are needed to better understand these sources of variability in the GSDs – for example, with a more systematic sampling of very fresh landslides where fines should not have been washed out. Thus, we encourage that such issues are anticipated in future studies and that field work is performed in a way which allows the spatial variability to be recorded. This would also enable future studies to include various landslide GSDs based on different assumptions or corrections. In this sense, collecting more measurements of landslides with both internal and superficial GSDs seems essential, especially when comparing young landslides with similar characteristics (e.g., lithology and height drop).

Figure 7Fraction of the landslide GSD that could be transported as both bed load and suspended load (a) and only as suspended load (b) as a function of the strength-normalized landslide drop height (as in Fig. 4) and river shear stress during a 10-year return flood (illustrated using different colors). Note that the suspended fractions at 300 and 380 Pa are identical in panel (b).


4.3 Implications for sediment transport in Taiwan

The landslide GSDs that we report contain mainly gravel, but they also contain a substantial fraction of boulders, which suggests that, after reaching floodplains and channels, the transport and evacuation of the material will require large floods. To compare these GSDs to typical shear stresses occurring in Taiwanese rivers, we use the shear stress map derived by Yanites et al. (2010b) from detailed measurement of the width, discharge and slope along the Peikang River. For a 10-year return flood with a discharge of 1000 m3 s−1, they found that shear stress, τ, mostly ranged from about 60 to about 380 Pa. For mountain channels with a typical gradient of about 2 % (Yanites et al.2010b), these shear stresses correspond to flood heights of between 0.3 and 2 m. To assess a threshold for bed load transport, we computed the grain size D for which the Shields number τ/(ρ/ρf-1)gD, where ρ and ρf are the grain and fluid density, respectively, was above a transport threshold of 0.045 (Lamb et al.2008). Similarly, for suspended load transport, we assessed the D value at which the shear velocity, estimated as U*=τ/ρf, was larger than the settling velocity of the grain Us as defined and calibrated by Ferguson and Church (2004) (Figs. 27). Even for an above average 10-year return flood, less than 25 % of most landslide deposits could be transported in suspension, except for LS-10 and LS-16 which had a suspended fraction of up to 50 %–70 %. When accounting for bed load transport, the largest shear stress of ∼380 Pa could not transport 5 %–25 % of the deposits for about half of the landslides, especially LS-6, LS-9n and LS-13. Considering smaller, but not uncommon, shear stresses (60–140 Pa) would result in an immobile fraction of 20 % to 40 % for most landslides, and up to 80 % for the three coarsest deposits.

Before discussing the implications of this, we highlight three main limitations that should be addressed by future work aiming at constraining the export of landslides deposits. First, the shear stress could not be adjusted to the local channel conditions in which the landslide occurred, neglecting specific width, discharge and gradient as well as relations between gradient and the critical Shields value (e.g., Lamb et al.2008) or the influence of landsliding on the channel itself (e.g., Kuo and Brierley2014). Second, we ignored armoring effects, in which a superficial layer of coarse grains inhibits the mobility of finer grains (Parker and Sutherland1990). In our case, considering armoring could particularly reduce transport for deposits where coarse grains are segregated at the toe or surface of the deposit, such as for LS-2, LS-3 or LS-8. Third, we did not consider debris flows and hyperconcentrated flows, which are frequent in Taiwan and sometimes reach the ocean (Dadson et al.2005; Lin et al.2005; Hsu et al.2010) and which would enhance sediment transport given their higher fluid density. Despite these sources of uncertainty, our results suggest that rapid evacuation of the sediment by suspension affects at most 30 % of most of the deposits in the relatively strong metasedimentary units, and most of the transport occurs as bed load. Further, only the largest (10-year return or more) floods will transport substantial parts of the deposit, meaning that large landslide events may load channels with a pulse of coarse sediments that require several decades to be evacuated. This is much longer than the transient pulse of enhanced landsliding (Marc et al.2015) and suspended sediment transport (Hovius et al.2011) observed after the Chi-Chi earthquake, which lasted less than 10 years, but is consistent with the ∼50-year timescales for enhanced lake sediment deposition (including bed load) after earthquakes observed in New Zealand (Howarth et al.2012; Wang et al.2020). This multidecadal timescale for sediment export seems consistent with the very large alluviation of river channels in southern Taiwan following intense flooding and landsliding triggered by Typhoon Morakot (Yanites et al.2018), which was still visible in 2015 (e.g., Taimali River) and at the time of writing in satellite imagery. Substantial aggradation, suspected to be long term, was also observed after the Chi-Chi earthquake (Yanites et al.2010a; Chen2009). More detailed modeling of the evacuation of landslide sediment (e.g., Yanites et al.2010a; Croissant et al.2017) could be combined with scenarios based on the detailed GSDs reported in this study to better quantify the dynamics and timescales of coarse sediment export after large landslides events.

5 Conclusions

We presented grain size distributions obtained from 17 landslide deposits in Taiwan. They have D50 and D84 values consistent with landslides reported in previous literature – between 15 and 200 mm and between 60 and 600 mm, respectively. We found that many deposits had significant spatial segregation in the downslope direction, with the lowest part of the deposits having 2–10 times coarser GSDs than the upper part of the deposits. For the three landslides in which we could sample the inner part of deposits, we also found that GSDs were 3–10 times finer than their surface counterparts. The presence and intensity of this segregation cannot be attributed to a single process, but kinetic sieving and deposit reworking are likely to play important roles. This internal variability could bias results obtained from local sampling of GSDs, such as sieve samples from a single pit.

Investigating the controls on landslide GSD variability, we observed a strong anticorrelation between the landslide drop height, width and inferred scar depth as well as the GSD percentiles for all of the landslides. Finer GSDs in the two landslides in unmetamorphosed, young sedimentary rocks can be well explained by normalizing the drop height by the rock strength. Further, modeling the source material median grain size with an exponential fining towards the surface, consistent with physical weathering models, we found that the reduction ratio from source material to landslide deposits matches the scaling proposed by Locat et al. (2006) and calibrated for rock avalanches with much larger volume and much higher point load strength than the ones we studied. Although future measurements on source rocks are needed for a complete demonstration, especially in terms of bedrock strength and fracture spacing, we suggest that simple geomorphic models coupling this fragmentation scaling with a model for regolith grain size (see Eq. 2) could provide a physically based first-order model for the GSD input to rivers by landslides in active orogens. Such an approach could be implemented into landscape evolution models accounting for sediment transport. Indeed, from our deposits, we also noted that even a 10-year flood may not be able to transport the coarsest fraction of many deposits, suggesting that floodplains and channels will likely need several decades to recover after large landslide events.

Data availability

The 28 GSDs (for the 17 landslides and each of their subsector samples) are available in the HydroShare open repository along with a shapefile comprising landslide locations and polygons derived from Google Earth: (Marc et al.2021).


The supplement related to this article is available online at:

Author contributions

OM designed the study and the field mission, performed all analyses, and wrote the paper. JMT and PM provided input for the field methodology and the result interpretations, and edited the paper. All authors participated in collecting the grain size data in the field.

Competing interests

The authors declare that they have no conflict of interest.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The grain count during the 2015 field campaign would not have been possible without the additional participation of Antonius Golly, Arnaud Burtin, Anne Schöpa and Niels Hovius, and they are warmly thanked for this. We are also grateful to Kristen Cook and Anne Schöpa for contributing pictures and to Sebastien Carretier for pointing us to the literature on physical weathering models. We wish to acknowledge Mikael Attal, one anonymous reviewer and the associate editor, Rebecca Hodge, for their constructive reviews that helped clarify and improve this paper. We also acknowledge Camille Noûs (last access: 13 August 2021), who embodies the broader scientific community and the collaborative and open nature of the creation and dissemination of knowledge, for its contribution to the design, methodology and interpretation of this study.

Financial support

The article processing charges for this open-access publication were covered by the Helmholtz Centre Potsdam – GFZ German Research Centre for Geosciences.

Review statement

This paper was edited by Rebecca Hodge and reviewed by Mikaël Attal and one anonymous referee.


Allen, P. A., Armitage, J. J., Whittaker, A. C., Michael, N. A., Roda-Boluda, D., and D'Arcy, M.: Fragmentation Model of the Grain Size Mix of Sediment Supplied to Basins, J. Geol., 123, 405–427,, 2015. a, b

Anderson, R. S., Anderson, S. P., and Tucker, G. E.: Rock damage and regolith transport by frost: an example of climate modulation of the geomorphology of the critical zone, Earth Surf. Proc. Land., 38, 299–316,, 2013. a, b, c, d, e

Armitage, J. J., Duller, R. A., Whittaker, A. C., and Allen, P. A.: Transformation of tectonic and climatic signals from source to sedimentary archive, Nat. Geosci., 4, 231–235,, 2011. a

Attal, M. and Lavé, J.: Changes of bedload characteristics along the Marsyandi River (central Nepal): Implications for understanding hillslope sediment supply, sediment load evolution along fluvial networks, and denudation in active orogenic belts, Geol. Soc. Am. Spec. Pap., 398, 143–171,, 2006. a, b, c, d, e, f

Attal, M., Mudd, S. M., Hurst, M. D., Weinman, B., Yoo, K., and Naylor, M.: Impact of change in erosion rate and landscape steepness on hillslope and fluvial sediments grain size in the Feather River basin (Sierra Nevada, California), Earth Surf. Dynam., 3, 201–222,, 2015. a, b

Benda, L. and Dunne, T.: Stochastic forcing of sediment routing and storage in channel networks, Water Resour. Res., 33, 2865–2880,, 1997. a

Campforts, B., Shobe, C. M., Steer, P., Vanmaercke, M., Lague, D., and Braun, J.: HyLands 1.0: a hybrid landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution, Geosci. Model Dev., 13, 3863–3886,, 2020. a

Carretier, S., Martinod, P., Reich, M., and Godderis, Y.: Modelling sediment clasts transport during landscape evolution, Earth Surf. Dynam., 4, 237–251,, 2016. a

Casagli, N., Ermini, L., and Rosati, G.: Determining grain size distribution of the material composing landslide dams in the Northern Apennines: sampling and processing methods, Eng. Geol., 69, 83–97,, 2003. a, b, c

Chau, K. T. and Wong, R. H. C.: Uniaxial compressive strength and point load strength of rocks, Int. J. Rock Mech. Min. Sci. Geomech. Abstr., 33, 183–188,, 1996. a

Chen, C.-Y.: Sedimentary impacts from landslides in the Tachia River Basin, Taiwan, Geomorphology, 105, 355–365,, 2009. a

Chung, C.-H. and Chang, F.-J.: A refined automated grain sizing method for estimating river-bed grain size distribution of digital images, J. Hydrol., 486, 224–233,, 2013. a, b

Clarke, B. A. and Burbank, D. W.: Quantifying bedrock-fracture patterns within the shallow subsurface: Implications for rock mass strength, bedrock landslides, and erodibility, J. Geophys. Res.-Earth, 116, F04009,, 2011. a, b

Cohen, S., Willgoose, G., and Hancock, G.: The mARM3D spatially distributed soil evolution model: Three-dimensional model framework and analysis of hillslope and landform responses, J. Geophys. Res.-Earth, 115, F04013,, 2010. a, b, c, d, e, f, g

Cook, K. L., Turowski, J. M., and Hovius, N.: A demonstration of the importance of bedload transport for fluvial bedrock erosion and knickpoint propagation, Earth Surf. Proc. Land., 38, 683–695,, 2013. a

Cook, K. L., Turowski, J. M., and Hovius, N.: River gorge eradication by downstream sweep erosion, Nat. Geosci., 7, 682–686,, 2014. a

Croissant, T., Lague, D., Steer, P., and Davy, P.: Rapid post-seismic landslide evacuation boosted by dynamic river width, Nat. Geosci., 10, 680–684,, 2017. a

Crosta, G. B., Frattini, P., and Fusi, N.: Fragmentation in the Val Pola rock avalanche, Italian Alps, J. Geophys. Res.-Earth, 112, F01006,, 2007. a, b, c, d, e

Dadson, S., Hovius, N., Pegg, S., Dade, W. B., Horng, M. J., and Chen, H.: Hyperpycnal river flows from an active mountain belt, J. Geophys. Res.-Earth, 110, F04016,, 2005. a

De Blasio, F. V. and Crosta, G. B.: Simple physical model for the fragmentation of rock avalanches, Acta Mechan., 225, 243–252,, 2014. a

Domej, G., Bourdeau, C., and Lenti, L.: Mean Landslide Geometries Inferred from a Global Database of Earthquake- and Non-Earthquake-Triggere Landslides, Ital. J. Eng. Geol. Environ., 17, 87–107,, 2017. a

Egholm, D. L., Knudsen, M. F., and Sandiford, M.: Lifespan of mountain ranges scaled by feedbacks between landsliding and erosion by rivers, Nature, 498, 475–478,, 2013. a

Emberson, R., Hovius, N., Galy, A., and Marc, O.: Chemical weathering in active mountain belts controlled by stochastic bedrock landsliding, Nat. Geosci., 9, 42–45,, 2016a. a

Emberson, R., Hovius, N., Galy, A., and Marc, O.: Oxidation of sulfides and rapid weathering in recent landslides, Earth Surf. Dynam., 4, 727–742,, 2016b. a

Ferguson, R. I. and Church, M.: A Simple Universal Equation for Grain Settling Velocity, J. Sediment. Res., 74, 933–937,, 2004. a

Gray, J. M. N. T.: Particle Segregation in Dense Granular Flows, Annu. Rev. Fluid Mech., 50, 407–433,, 2018. a

Guerit, L., Barrier, L., Narteau, C., Métivier, F., Liu, Y., Lajeunesse, E., Gayer, E., Meunier, P., Malverti, L., and Ye, B.: The Grain-size Patchiness of Braided Gravel-Bed Streams – example of the Urumqi River (northeast Tian Shan, China), Adv. Geosci., 37, 27–39,, 2014. a

Guerit, L., Barrier, L., Liu, Y., Narteau, C., Lajeunesse, E., Gayer, E., and Métivier, F.: Uniform grain-size distribution in the active layer of a shallow, gravel-bedded, braided river (the Urumqi River, China) and implications for paleo-hydrology, Earth Surf. Dynam., 6, 1011–1021,, 2018. a

Hovius, N., Stark, C. P., and Allen, P. A.: Sediment flux from a mountain belt derived by landslide mapping, Geology, 25, 231–234,<0231:SFFAMB>2.3.CO;2, 1997. a

Hovius, N., Stark, C. P., Chu, H., and Lin, J.: Supply and Removal of Sediment in a Landslide Dominated Mountain Belt: Central Range, Taiwan, J. Geol., 108, 73–89,, 2000. a

Hovius, N., Meunier, P., Lin, C., Chen, H., Chen, Y., Dadson, S., Horng, M., and Lines, M.: Prolonged seismically induced erosion and the mass balance of a large earthquake, Earth Planet. Sc. Lett., 304, 347–355,, 2011. a

Howarth, J. D., Fitzsimons, S. J., Norris, R. J., and Jacobsen, G. E.: Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault, New Zealand, Geology, 40, 1091–1094,, 2012. a

Hsu, S. M., Chiou, L. B., Lin, G. F., Chao, C. H., Wen, H. Y., and Ku, C. Y.: Applications of simulation technique on debris-flow hazard zone delineation: a case study in Hualien County, Taiwan, Nat. Hazards Earth Syst. Sci., 10, 535–545,, 2010. a

Ibbeken, H.: Jointed source rock and fluvial gravels controlled by Rosin's law; a grain-size study in Calabria, South Italy, J. Sediment. Res., 53, 1213–1231,, 1983. a, b

Keefer, D. K.: Landslides caused by earthquakes, Geol. Soc. Am. Bull., 95, 406–421,<406:LCBE>2.0.CO;2, 1984. a

Kellerhals, R. and Bray, D. I.: Sampling Procedures for Coarse Fluvial Sediments, J. Hydraul. Div., 97, 1165–1180, 1971. a

Kuo, C.-W. and Brierley, G.: The influence of landscape connectivity and landslide dynamics upon channel adjustments and sediment flux in the Liwu Basin, Taiwan, Earth Surf. Proc. Land., 39, 2038–2055,, 2014. a

Lamb, M. P., Dietrich, W. E., and Venditti, J. G.: Is the critical Shields stress for incipient sediment motion dependent on channel-bed slope?, J. Geophys. Res.-Earth, 113, F02008,, 2008. a, b

Larsen, I., Montgomery, D., and Korup, O.: Landslide erosion controlled by hillslope material, Nat. Geosci., 3, 247–251, 2010. a

Lin, C.-P., Wang, Y.-M., Tfwala, S. S., and Chen, C.-N.: The Variation of Riverbed Material due to Tropical Storms in Shi-Wen River, Taiwan, Scient. World J., 2014, e580936,, 2014. a

Lin, G.-W., Chen, H., Hovius, N., Horng, M.-J., Dadson, S., Meunier, P., and Lines, M.: Effects of earthquake and cyclone sequencing on landsliding and fluvial sediment transfer in a mountain catchment, Earth Surf. Proc. Land., 33, 1354–1373,, 2008. a, b, c, d

Lin, M.-L., Wang, K.-L., and Huang, J.-J.: Debris flow run off simulation and verification - case study of Chen-You-Lan Watershed, Taiwan, Nat. Hazards Earth Syst. Sci., 5, 439–445,, 2005. a

Locat, P., Couture, R., Leroueil, S., Locat, J., and Jaboyedoff, M.: Fragmentation energy in rock avalanches, Can. Geotech. J., 43, 830–851,, 2006. a, b, c, d, e, f, g, h, i, j

Marc, O., Hovius, N., Meunier, P., Uchida, T., and Hayashi, S.: Transient changes of landslide rates after earthquakes, Geology, 43, 883–886,, 2015. a

Marc, O., Stumpf, A., Malet, J.-P., Gosset, M., Uchida, T., and Chiang, S.-H.: Initial insights from a global database of rainfall-induced landslide inventories: the weak influence of slope and strong influence of total storm rainfall, Earth Surf. Dynam., 6, 903–922,, 2018. a, b

Marc, O., Behling, R., Andermann, C., Turowski, J. M., Illien, L., Roessner, S., and Hovius, N.: Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides, Earth Surf. Dynam., 7, 107–128,, 2019. a, b, c

Marc, O., Turowski, J. M., and Meunier, P.: Grain Size Distribution of 17 Taiwanese landslide deposits, [data set], Hydroshare,, 2021. a

Marshall, J. A. and Sklar, L. S.: Mining soil databases for landscape-scale patterns in the abundance and size distribution of hillslope rock fragments, Earth Surf. Proc. Land., 37, 287–300,, 2012. a

Meunier, P., Hovius, N., and Haines, J. A.: Topographic site effects and the location of earthquake induced landslides, Earth Planet. Sc. Lett., 275, 221–232,, 2008. a

Neely, A. B. and DiBiase, R. A.: Drainage Area, Bedrock Fracture Spacing, and Weathering Controls on Landscape-Scale Patterns in Surface Sediment Grain Size, J. Geophys. Res.-Earth, 125, e2020JF005560,, 2020. a, b

Neely, A. B., DiBiase, R. A., Corbett, L. B., Bierman, P. R., and Caffee, M. W.: Bedrock fracture density controls on hillslope erodibility in steep, rocky landscapes with patchy soil cover, southern California, USA, Earth Planet. Sc. Lett., 522, 186–197,, 2019. a

Nishiguchi, Y., Uchida, T., Takezawa, N., Ishizuka, T., and Mizuyama, T.: Runout Characteristics and Grain Size Distribution of Large-scale Debris Flows Triggered by Deep Catastrophic Landslides, Int. J. Erosion Control Eng., 5, 16–26,, 2012. a

Parker, G. and Sutherland, A. J.: Fluvial armor, J. Hydraul. Res., 28, 529–544,, 1990. a

Rault, C., Robert, A., Marc, O., Hovius, N., and Meunier, P.: Seismic and geologic controls on spatial clustering of landslides in three large earthquakes, Earth Surf. Dynam., 7, 829–839,, 2019. a

Riebe, C. S., Sklar, L. S., Lukens, C. E., and Shuster, D. L.: Climate and topography control the size and flux of sediment produced on steep mountain slopes, P. Natl. Acad. Sci. USA, 112, 15574–15579,, 2015. a

Roda‐Boluda, D. C., D'Arcy, M., McDonald, J., and Whittaker, A. C.: Lithological controls on hillslope sediment supply: insights from landslide activity and grain size distributions, Earth Surf. Proc. Land., 43, 956–977,, 2018. a, b, c

Ruiz-Carulla, R. and Corominas, J.: Analysis of Rockfalls by Means of a Fractal Fragmentation Model, Rock Mech. Rock Eng., 53, 1433–1455,, 2020. a, b

Ruiz-Carulla, R., Corominas, J., and Mavrouli, O.: A methodology to obtain the block size distribution of fragmental rockfall deposits, Landslides, 12, 815–825,, 2015. a, b, c, d

Savage, S. B. and Lun, C. K. K.: Particle size segregation in inclined chute flow of dry cohesionless granular solids, J. Fluid Mech., 189, 311–335,, 1988. a

Sklar, L. S. and Dietrich, W. E.: Sediment and rock strength controls on river incision into bedrock, Geology, 29, 1087–1090,<1087:SARSCO>2.0.CO;2, 2001. a

Sklar, L. S., Riebe, C. S., Marshall, J. A., Genetti, J., Leclere, S., Lukens, C. L., and Merces, V.: The problem of predicting the size distribution of sediment supplied by hillslopes to rivers, Geomorphology, 277, 31–49,, 2017. a, b, c, d

Stephens, M. A.: EDF Statistics for Goodness of Fit and Some Comparisons, J. Am. Stat. Assoc., 69, 730–737,, 1974. a

Struck, M., Andermann, C., Hovius, N., Korup, O., Turowski, J. M., Bista, R., Pandit, H. P., and Dahal, R. K.: Monsoonal hillslope processes determine grain size-specific suspended sediment fluxes in a trans-Himalayan river: Mass wasting determines sediment caliber, Geophys. Res. Lett., 42, 2302–2308,, 2015. a

Tanyaş, H., van Westen, C. J., Allstadt, K. E., Anna Nowicki Jessee, M., Görüm, T., Jibson, R. W., Godt, J. W., Sato, H. P., Schmitt, R. G., Marc, O., and Hovius, N.: Presentation and Analysis of a Worldwide Database of Earthquake-Induced Landslide Inventories, J. Geophys. Res.-Earth, 122, 2017JF004236,, 2017. a

Turowski, J. M.: Alluvial cover controlling the width, slope and sinuosity of bedrock channels, Earth Sur. Dynam., 6, 29–48,, 2018. a

Vallance, J. W. and Savage, S. B.: Particle Segregation in Granular Flows Down Chutes, in: IUTAM Symposium on Segregation in Granular Flows, edited by: Rosato, A. D. and Blackmore, D. L., Springer Netherlands, Dordrecht, 31–51, 2000. a

Varnes, D.: Slope movement types and processes, in: Landslides, analysis and control, no. 170 in Special Report Transportation Research Board, edited by: Schuster, R. L. and Krizek, R. J., National Academy of Science, Washington, 11–33, 1978. a

Verdian, J. P., Sklar, L. S., Riebe, C. S., and Moore, J. R.: Sediment size on talus slopes correlates with fracture spacing on bedrock cliffs: Implications for predicting initial sediment size distributions on hillslopes, Earth Surf. Dynam. Discuss. [preprint],, in review, 2020.  a

Wang, J., Howarth, J. D., McClymont, E. L., Densmore, A. L., Fitzsimons, S. J., Croissant, T., Gröcke, D. R., West, M. D., Harvey, E. L., Frith, N. V., Garnett, M. H., and Hilton, R. G.: Long-term patterns of hillslope erosion by earthquake-induced landslides shape mountain landscapes, Sci. Adv., 6, eaaz6446,, 2020. a

Whittaker, A. C., Duller, R. A., Springett, J., Smithells, R. A., Whitchurch, A. L., and Allen, P. A.: Decoding downstream trends in stratigraphic grain size as a function of tectonic subsidence and sediment supply, GSA Bull., 123, 1363–1382,, publisher: GeoScienceWorld, 2011. a

Yanites, B., Tucker, G., Mueller, K., and Chen, Y.: How rivers react to large earthquakes: Evidence from central Taiwan, Geology, 38, 639–642, 2010a. a, b

Yanites, B. J., Tucker, G. E., Mueller, K. J., Chen, Y.-G., Wilcox, T., Huang, S.-Y., and Shi, K.-W.: Incision and channel morphology across active structures along the Peikang River, central Taiwan: Implications for the importance of channel width, GSA Bull., 122, 1192–1208,, 2010b. a, b

Yanites, B. J., Mitchell, N. A., Bregy, J. C., Carlson, G. A., Cataldo, K., Holahan, M., Johnston, G. H., Nelson, A., Valenza, J., and Wanker, M.: Landslides control the spatial and temporal variation of channel width in southern Taiwan: Implications for landscape evolution and cascading hazards in steep, tectonically active landscapes, Earth Surf. Proc. Land., 43, 1782–1797,, 2018. a

Short summary
The size of grains delivered to rivers is an essential parameter for understanding erosion and sediment transport and their related hazards. In mountains, landslides deliver these rock fragments, but few studies have analyzed the landslide properties that control the resulting sizes. We present measurements on 17 landslides from Taiwan and show that their grain sizes depend on rock strength, landslide depth and drop height, thereby validating and updating a previous theory on fragmentation.