Denudation Rate Changes Along a Fast‐Eroding Mountainous River With Slate Headwaters in Taiwan From 10Be (Meteoric)/9Be Ratios

The novel 10Be (meteoric)/9Be system, where 10Be is delivered by precipitation and stable 9Be is released by weathering, provides denudation rates over weathering‐erosion timescales. The new tool is applicable to quartz‐poor lithologies, for example, mafic rock and claystone, which are not readily accessible by the commonly used in situ‐produced 10Be in quartz. We provide a first application of this proxy to a tectonically active mountainous river, the Zhuoshui River in Taiwan. Taiwan Rivers supply a disproportionately high suspended and dissolved flux to the oceans and are often underlain by fine‐grained shale/slate. 10Be (meteoric)/9Be‐derived denudation rates (Dmet) from the Zhuoshui Catchment are highest in the slate‐dominated headwaters (4–8 mm/year), and much lower (1–2 mm/year) along the midlower reaches with mixed lithologies. At the basin‐wide scale, we find a poor correlation between Dmet and basin‐averaged channel steepness despite a small climatic gradient. Because large lithological heterogeneities exist in this basin, we invoke a lithological effect to explain this poor correlation. Relying on a revised stream power incision model that incorporates rock erodibility, the resulting lithology‐ and runoff‐adjusted ksn (kLrsn) can be reconciled with denudation rates with the highest erodibility predicted to prevail in the Miocene slate of low metamorphic grade and high fracture density. This model suggests that the lithological heterogeneity can alter the coupling between surface denudation and channel morphology. On a broader perspective, the successful application of the 10Be (meteoric)/9Be proxy shows its applicability as a tracer for erosion and sediment transport processes in fast‐eroding mountain belts underlain by slate lithologies.


Introduction
Small mountainous rivers (area <10 4 km 2 and maximum elevation >1 km) discharging to the oceans have been in the focus of sediment flux studies for many years because of their disproportionately high sediment discharge (Hovius et al., 1997;Milliman & Syvitski, 1992). As a fast-eroding mountain belt, the Taiwan orogen alone bears 8 of 13 rivers in the world with sediment yields higher than 10 4 t/km 2 /year (Milliman & Farnsworth, 2011). The extremely high erosion rates of Taiwan catchments are mainly attributed to active tectonics and frequent typhoon attacks (Hovius & Lin, 2000).
To invoke controls on modern topography and landscape evolution of mountain belts, millennial-scale denudation rates are commonly determined using in situ-produced cosmogenic nuclides (e.g., Bierman et al., 2002;Granger & Schaller, 2014;von Blanckenburg, 2005). In most cases, in situ 10 Be in quartz is the first choice as production rates and scaling are well known (Balco et al., 2008). In Taiwan, basin-wide in situ 10 Be-derived denudation rates have been obtained for a number of rivers (Derrieux et al., 2014;Fellin et al., 2017;Siame et al., 2011). From these studies, denudation rates of 1-3 mm/year have been estimated for the prowedge of the mountain belt (West Taiwan), and higher denudation rates on the order of 4-5 mm/year for the retro-wedge (East Taiwan; Derrieux et al., 2014). However, considering that quartz-poor slate dominates the alpine region of Taiwan (Garzanti & Resentini, 2016), in situ 10 Be-derived denudation rates may not fully represent basin-wide erosion processes. In this case, a recently developed denudation proxy, the 10 Be (meteoric)/ 9 Be ratio, provides an alternative choice, as it does not rely on quartz lithologies.
Meteoric 10 Be is produced at an approximately known flux in the atmosphere and accumulates at the Earth's surface (Willenbring & von Blanckenburg, 2010), while the release rate of stable 9 Be from bedrocks depends on chemical weathering. Von  developed a series of steady-state mass balance equations that describe the 10 Be (meteoric)/ 9 Be system from a soil profile up to the whole catchment. Compared to the sister nuclide in situ 10 Be, the advantages of meteoric 10 Be lie in its higher concentrations and thus smaller sample amounts needed (~1 g) for measurement, and its applicability to quartz-free lithologies . Although mineralogical sorting, sorption capacity, and retentivity affect the distribution of meteoric 10 Be and 9 Be concentrations (Singleton et al., 2016;You et al., 1989), the 10 Be (meteoric)/ 9 Be ratio is largely independent of these processes  when Be is chemically extracted from the reactive phase . Thus, the 10 Be (meteoric)/ 9 Be ratio has been shown to record basin-wide denudation rates in the Amazon basin (Wittmann et al., 2015), the Ganga basin (Rahaman et al., 2017), the Potomac Basin (Portenga et al., 2019), and very small catchments (<1 km 2 ) draining mafic and ultramafic rocks in the Slavkov Forest (Dannhaus et al., 2018). Even though in these studies denudation rates derived from in situ 10 Be and 10 Be (meteoric)/ 9 Be ratios mostly agree within a factor of 2, larger differences between the two methods remain that might be due to methodological causes or reflect true variability in geologic conditions to which both methods respond differently. Hence, further comparison between both approaches is needed to show the applicability of 10 Be (meteoric)/ 9 Be ratios across landscapes. Specifically, the 10 Be (meteoric)/ 9 Be system has not been applied to small tectonically active mountainous rivers that contribute a large proportion of global sediment and dissolved fluxes (Milliman & Farnsworth, 2011).
In this study, we focus on the longest river in Taiwan, that is, the Zhuoshui (Choshui) River, which drains slate lithologies in its headwaters and sandstone-dominated lithologies along its lower reaches. This river with high sediment yield (>10 4 t/km 2 /year) can deliver the decadal-averaged annual sediment load during a single typhoon event (Goldsmith et al., 2008;Milliman et al., 2017). Sediment samples were collected from the mountain headwaters to the outlet in two different years ( Figure 1). We compare our new 10 Be (meteoric)/ 9 Be-derived denudation rates with published denudation and erosion data in Taiwan derived from other methods, including sediment gauging, in situ 10 Be, and detrital thermochronology. Our aim is to investigate the variations in denudation from the headwaters to the river mouth during sediment transport and mixing, and the dependence of denudation on topographic and lithological factors. This study will provide significant insights on using cosmogenic nuclide tools to study long-term denudation processes in fast-eroding mountain belts with exposure of fine-grained, quartz-poor lithologies, including some small mountainous catchments in high-standing islands of Southeastern Asia and Oceania (Jiao et al., 2017;Liu et al., 2016).

Topographic State and Denudation of the Taiwan Orogen
The Taiwan orogen was established by the oblique collision between the Luzon Arc and the Asian continental margin since~6 Ma, with the range gradually propagating southward (Huang et al., 2006). Given the southward propagation of mountain building processes, the temporal evolution of the topographic state in the orogen can be detected from the latitudinal pattern of multiple topographic metrics, for example, mean elevation, local slope, and channel steepness (Chen & Willett, 2016;Stolar et al., 2007). In general, the values of these metrics are relatively stable and high in the central part (e.g., 125-275 km from the southern tip of this island), indicating the mature stage of arc-continent collision and plausible topographic steady state, that is, the balance between long-term uplift and denudation. In comparison, these metrics gradually decrease southward and northward of the mature collision zone, corresponding to initial and postcollision stages (Chen & Willett, 2016;Stolar et al., 2007).
Denudation of the Taiwan orogen has been investigated in detail with various methods spanning different timescales (see synthesis in Dadson et al., 2003, andFellin et al., 2017). Short-term denudation rates have been quantified by sediment gauging over single events (e.g., one typhoon event; Goldsmith et al., 2008;Milliman et al., 2017) or over decadal scale (Kao & Milliman, 2008) and can be apportioned to diverse sediment source rocks when coupled with petrographic and mineralogical fingerprints (Resentini et al., 2017). Longer-term methods include reservoir sedimentary records (Montgomery et al., 2014) over centennial scales, in situ-produced cosmogenic 10 Be measured in quartz of detrital sediments over the timescale of several hundreds to thousands of years (Derrieux et al., 2014;Fellin et al., 2017;Siame et al., 2011), and also low-  Table S1.

10.1029/2019JF005251
Journal of Geophysical Research: Earth Surface temperature thermochronometry over millions of years (Hsu et al., 2016;Lee et al., 2006). The integration timescales of these methods relate to different controlling factors of denudation. Over short timescales, denudation rates are largely influenced by the frequency of typhoons and earthquakes Hovius et al., 2011). Over long timescales, plate motion and tectonic deformation velocities dominate denudation rates (Hsu et al., 2016). However, even though denudation rates vary greatly for each method and integration timescale, they are overall within the same order of 1-10 mm/year.

Geology and Hydrology
The Taiwan orogen is mainly composed of (meta-)sedimentary rocks ranging in age from the Paleozoic to the Quaternary. Six geological units can be identified from west to east: Coastal Plain, Western Foothills, Hsuehshan Range, Backbone Range, Tananao Complex (Tailuko Belt and Yuli Belt), and Coastal Range. The Backbone Range and the Tananao Complex can be grouped together as the Central Range ( Figure 1). Regarding exposed strata in the Zhuoshui Basin (Figure 1), the Miocene and Eocene slate (the Backbone Range) and the Eocene-upper Oligocene slate/metasandstone (the Hsuehshan Range) outcrop in the upper and middle reaches, and the Miocene-Pleistocene sandstone/shale (the Western Foothills) outcrops in the lower reaches. The narrow Coastal Plain in the downstream part of this basin is characterized by Quaternary alluvial deposits (Figure 1), comprising~4% of the total basin area. Detailed lithological information of each tectonic unit is provided in Table S1 in the supporting information.
As the longest river in Taiwan, the Zhuoshui River has a total length of~200 km and a basin area of 3×10 3 km 2 . It originates from the Central Range and flows westward to the Taiwan Strait, with the contributions of two major tributaries, the Chenyoulan River (area of~450 km 2 ) and the Qingshui River (area of 410 km 2 ; Figure 1). According to sediment gauging records covering several decades, the decadal-scale mean suspended sediment load of the Zhuoshui River is estimated at 40 Mt/year (Kao & Milliman, 2008) or 54 Mt/year (Dadson et al., 2003) based on different processing methods. However, sediment gauging observations conducted during several typhoon events show that 61-70 Mt of suspended load can also be transported within several days (Goldsmith et al., 2008;Milliman et al., 2017), indicating extremely rapid mass wasting.

Sampling Details
We conducted two sampling campaigns in nonflood seasons along the Zhuoshui River ( Figure 1a). Seven bedload samples (mostly sand) were collected in May 2013, including two tributary samples from the Chenyoulan River and the Qingshui River. In early June 2017, six bedload samples (mostly sand) were taken along the main stem (Table 1). All these samples were collected near the water line to record recent transport events. The upstream basins of all the sampling points are larger than 100 km 2 , to average-out potential bias introduced on cosmogenic 10 Be by unrepresentative distribution of recent landslides (Yanites et al., 2009). In addition, one delta-facies sample with modern age was taken from the top (depth:~1.1 m) of a~104-m-long sediment core (JRD-S) at the Zhuoshui River mouth (Yang et al., 2017). The purpose of collecting sediment samples between different years is to test the consistency in the downstream trend of Be isotopes, which may be of importance in such a landscape frequently influenced by stochastic events (Deng et al., 2019). Note that this data set does not aim at strictly testing the temporal replicability, which requires temporal replicates at each location (Gonzalez et al., 2017). Furthermore, we collected 14 bedrock samples from different lithological units contributing to Zhuoshui River sediment (see Table S2), for characterizing the 9 Be concentration of every lithological unit.
3.2. Calculation of Denudation Rates Using 10 Be (Meteoric)/ 9 Be Ratios 3.2.1. Conceptual Framework In order to quantify Earth surface erosion and weathering, a steady-state mass balance for the 10 Be (meteoric)/ 9 Be system was developed by von . One primary assumption for this framework is that the inputs and outputs of both isotopes are balanced in the drainage basin over weathering and erosion timescales. For the 10  (2) J 10 Be atm includes the atmospheric depositional flux (F 10 Be met ; in atoms/m 2 /year) over a certain drainage area A riv,i (in m 2 ) for every subcatchment, and the 10 Be depositional flux is delivered to Earth surface by wet and dry deposition (Willenbring & von Blanckenburg, 2010). J 10 Be riv gives the flux of particulate-bound 10 Be ([ 10 Be] reac , in atoms/kg) eroded from the same drainage area at an independently known erosion rate (E, in kg/m 2 / year). The reactive 10 Be can be leached from sediment by sequential chemical extraction . This reactive Be is absorbed onto mineral surfaces or coprecipitated into secondary minerals, comprising mainly amorphous (called am-ox) and crystalline (called x-ox) Fe-Al-oxyhydroxides. At low pH, Be is soluble (see below) and dissolved 10 Be ([ 10 Be] diss , in atoms/L) in river water should be accounted for when assessing the output of 10 Be from a catchment. The ratio of J 10 Be riv to J 10 Be atm is similar to the erosion index developed by Brown et al. (1988), and a J 10 Be riv /J 10 Be atm ratio of 1 indicates the balance between input and output 10 Be fluxes. Regarding the stable 9 Be system, 9 Be contained in the unweathered bedrock (termed [ 9 Be] parent , in mg/kg) is partitioned by denudation (the sum of erosion and weathering) into three components: reactive 9 Be ([ 9 Be] reac , in mg/kg), dissolved 9 Be ([ 9 Be] diss , in mg/L), and residual 9 Be remaining in primary minerals ([ 9 Be] min , in mg/kg). In large catchments, where the lithologic composition approaches that of the average upper continental crust, [ 9 Be] parent is inferred to be 2.5 ± 0.5 ppm . In small (e.g., creek-scale, area <1 km 2 ) river basins, [ 9 Be] parent needs to be obtained from measuring local bedrock samples (Dannhaus et al., 2018). In this study, we use an area-weighted 9 Be concentration from all the geological units in the catchment to represent [ 9 Be] parent (see Table S2).
The dissolved Be ( 10 Be and 9 Be) flux is often not available. Instead, the pH-dependent Be partition coefficient K d ([Be] reac /[Be] diss , in L/kg) can be used to estimate the distribution between [Be] reac and The samples are sorted according to the upstream distance from the river mouth. "ZS13-" and "ZS17-" indicate samples collected in 2013 (in bold) and 2017, respectively. JRD-S1 is the core sediment at the river mouth. Samples with asterisk (*) are from tributaries, and ZS13-7 from the Chenyoulan River and ZS13-10 from the Qingshui River. b The normalized channel steepness index (k sn ) is calculated using a concavity index of 0.5 and a critical threshold drainage area (A cr ) of 1 km 2 following Fellin et al. (2017). The basin-averaged k sn is then calculated as mean k sn of all individual reaches (1-km segments). c The basinaveraged annual precipitation data set (with the observation period of 1970-2000) is sourced from WorldClim Version 2 (Fick & Hijmans, 2017). d Area-weighted 9 Be concentration of parent bedrock is calculated based on the average value and standard deviation of bedrock 9 Be data from each geological unit (see Table S2).

10.1029/2019JF005251
Journal of Geophysical Research: Earth Surface [Be] diss and hence estimate the potential loss of reactive Be due to desorption. In addition to pH values, the ratio of river water flux (or runoff R, in m/year) to sediment load (or erosion rate E, in kg/m 2 /year) affects the partitioning of the Be export flux into the dissolved phase. When R/E ≪ K d , the loss of reactive Be to the dissolved fraction is negligible and can be ignored . The pH range of the Zhuoshui River water measured in June 2017 is 7.5-9.4, which corresponds to a minimum K d value of~3.2×10 5 L/kg based on experimentally derived K d -pH relationships (You et al., 1989). In comparison, gauging data of water and suspended load from recent decades (Kao & Milliman, 2008) show that the average R/E reaches~89 L/kg and is hence almost 4 orders of magnitude lower than the K d value.
Since the partitioning of Be to the dissolved fraction is negligible in the Zhuoshui Basin, the simplified equation for calculating denudation rates (D met in m/year, using a sediment density ρ sedi of 2.65×10 3 kg/m 3 ), derived from the input-output balance of 10 Be and 9 Be, can be used: The detailed derivation for equation (3) is given in von . Different from the single isotope ( 10 Be) system, ( 10 Be/ 9 Be) reac in equation (3) is independent of mineralogical sorting, sorption capacity, and Be retentivity . Still, the term ([ 9 Be] min /[ 9 Be] reac ) in equation (3) may be influenced by the above processes to some extent, because reactive 9 Be is potentially enriched over [ 9 Be] min in fine sediments Wittmann et al., 2012). To minimize this bias, integration of river depth profiles is the most precise method . However, sampling river depth profiles may not be the most suitable strategy for such a dynamic mountainous river like the Zhuoshui River. Low water depth in dry seasons and highly turbulent flow during typhoon events may preclude the development of a well-sorted profile to integrate the whole grain-size range (Deng et al., 2019). Considering that silt-and clay-sized suspended sediment should dominate the total sediment load in the Zhuoshui River Turowski et al., 2010), analyzing a narrow grain-size fraction of bedload that is similar to the grain-size range of suspended load ensures a relatively uniform grain-size distribution and will largely provide the Be isotopic signal of the total sediment load.

Analytical Methods
After overnight oven-drying, river sediment samples were dry-sieved to the grain-size range of 30-63 μm, as one of the major fractions of suspended load (Milliman et al., 2017). Next, the extraction of 10 Be and 9 Be was performed based on the methods described in Wittmann et al. (2012). About 1 g of the sieved sediment was treated subsequently with 0.5 M HCl for extracting the amorphous oxide fraction (am-ox), and then with 1 M hydroxylamine-hydrochloride solution for the crystalline oxide fraction (x-ox). After these sequential extraction steps, the mineral residue (min) was decomposed with a mixture of hydrofluoric acid and aqua regia. The total reactive fraction was obtained by summing up both the am-ox and x-ox fractions.
We split each of these fractions (i.e., am-ox, x-ox, and min). One split was analyzed for stable 9 Be concentrations by Inductively Coupled Plasma-Optical Emission Spectroscopy (ICP-OES, Varian 720-ES). The second split of am-ox and x-ox fractions was spiked with a 9 Be carrier of known weight and purified for 10 Be analysis according to established methods (e.g., Wittmann et al., 2012). 10 Be concentrations were obtained from accelerator mass spectrometry measurements of 10 Be/ 9 Be ratios (with carrier) at University of Cologne, and the standard KN01-6-2 with a 10 Be/ 9 Be ratio of 5.35 × 10 -13 was applied (Dewald et al., 2013).
For river sediments collected in 2013 and the core sediment, the am-ox and x-ox fractions were analyzed separately; for river sediments sampled in 2017, am-ox and x-ox fractions were combined as reactive fractions before measurement and then analyzed for 10 Be and 9 Be. Essentially, there is no difference in Be isotopic data of the reactive fraction between two kinds of measurements.
Regarding bedrock samples, we used the chemical digestion procedure described by Deng et al. (2019). 9 Be concentrations of bulk samples were obtained by ICP-OES (Varian 720-ES) at the GFZ German Research Centre for Geosciences or ICP-MS (Agilent 7900) at the State Key Laboratory of Marine Geology, Tongji University.

Extraction of Topographic Metrics
For topographic analysis, we extracted the basin-averaged slope, local relief, and channel steepness index using the MATLAB-based software TopoToolbox 2 (Schwanghart & Scherler, 2014) based on SRTM (Shuttle Radar Topography Mission) 30-m resolution Digital Elevation Model (http://gdex.cr.usgs.gov/). Local slope for each pixel in the Digital Elevation Model is the steepest downward gradient according to neighborhood of eight connected pixels, and the basin-averaged slope is simply the average slope of all the pixels in the whole basin. Relief is a scale-dependent measurement, and local relief of each pixel is the elevation range within a specific radius (5 km in this case). Basin-averaged local relief is the average local relief of all the pixels within a catchment.
Furthermore, the channel steepness index is commonly used as a key metric characterizing channel morphology and derived from the power law correlation between local channel slope (S) and upstream drainage area (A; Flint, 1974] according to where k s and θ are referred to as the local channel steepness and concavity index, respectively. θ commonly falls in a narrow range (0.35-0.65; Wobus et al., 2006) and does not vary with rock uplift rate or lithology in a statistically significant way (Lague, 2014). This scaling relationship only holds for a drainage area larger than a critical threshold (A cr ), interpreted as a transition from colluvial channels to fluvial channels (Kirby & Whipple, 2012). To establish a link between channel incision and morphology in detachment-limited catchments, a commonly used approach is the stream power incision model (SPIM) as below (Whipple & Tucker, 1999): where D is the channel incision rate or denudation rate over the basin scale, K is the dimensional erosional efficiency related to lithology, climate and river morphology, and m and n are positive constants. θ in equation (4) can be taken as the m/n ratio if channels are in equilibrium and D is spatially uniform (Lague, 2014;Whipple & Tucker, 1999). Next, by combining equations (4) and (5), the local channel steepness can be described as a function of channel incision rate (Wobus et al., 2006): The positive correlation between k s and D holds when climatic factors and substrate rock properties (related to K) are relatively stable. Furthermore, at topographic steady state, the incision rate is equal to the rock uplift rate (U) by definition, and thus, k s is also considered as a proxy for uplift rate (Kirby & Whipple, 2012). In practice, θ and k s can be obtained from the linear regression between logs S and A (equation (4)), as slope and intercept, respectively (Wobus et al., 2006). However, since k s is very sensitive to the selection of θ, a fixed reference concavity index is commonly used, and a normalized local channel steepness index (k sn = SA θ ref Þ is calculated (Kirby & Whipple, 2012). As such, reach-scale k sn can be compared within a basin and linked to local variability in uplift rates or erosional efficiency (equation (6)). We adopt a θ ref of 0.5 in this study, as the concavity index of river channels in East Taiwan mainly varies between 0.4 and 0.6 . The distribution of normalized local channel steepness is shown in Figure 1b. The basin-averaged k sn is then calculated as mean k sn of all individual reaches (1-km segments).

Results
The spatial variations of [ 10 Be] reac , [ 9 Be] reac and ( 10 Be/ 9 Be) reac along the river are shown in Figure 2 and Table 2. [ 10 Be] reac generally increases downstream for samples collected in both 2013 and 2017, with a range of 3.6 × 10 5 −4.7 × 10 6 atoms/g. The average [ 10 Be] reac in the upper main stem (defined as upstream distance from the river mouth >110 km, 5.3 × 10 5 atoms/g, n = 4) is significantly lower than that along the floodplain reaches (upstream distance <48 km, 2.6 × 10 6 atoms/g, n = 5). The average [ 10 Be] reac of tributary samples from the Chenyoulan and Qingshui Rivers (3.9 × 10 6 atoms/g) is approximately 1.5-fold higher than that in the downstream floodplain. [ 9 Be] reac of sediment samples varies from 0.22 to 0.62 μg/g. In general, [ 9 Be] reac is slightly higher in upper reaches (0.45 μg/g on average, n = 4) than in the downstream floodplain (0.32 μg/g on average, n = 5). For samples where both am-ox and x-ox fractions were obtained (Table 2), the am-ox fraction dominates the reactive fraction for both 10 Be and 9 Be concentrations, accounting for 70-92%.
The bedrock 9 Be concentration in the upper reaches is higher than that in the lower reaches, ranging from 2.8 to 2.9 μg/g in the Backbone Range and Hsuehshan Range to 1.2 μg/g in the Western Foothills (see Table  S2). The [ 9 Be] parent (in equation (3)) of each sample (Table 1) can thus be calculated as an area-weighted 9 Be concentration of all the geological units. The resulting [ 9 Be] parent decreases from upper main stem (2.9 μg/g) to downstream floodplain (2.3 μg/g). The two tributaries, the Chenyoulan and Qingshui Rivers, contain relatively low [ 9 Be] parent of 2.0 and 1.2 μg/g, respectively.  Note.
(1) The samples are sorted according to the upstream distance from the river mouth.
(2) Measured chemical fractions include the amorphous oxide-bound fraction (am-ox), crystalline oxidebound fraction (x-ox), reactive fraction (reac, the sum of am-ox and x-ox), and silicate residue fraction (min). (3) (1)), the basin-wide meteoric 10 Be depositional flux F 10 Be met needs to be known. F 10 Be met can be obtained by several methods, including long-term meteoric 10 Be inventories measured in independently dated soils (Ouimet et al., 2015), short-term measurements of 10 Be in precipitation (Graly et al., 2011), or atmospheric general circulation models (GCMs) for global meteoric 10 Be deposition (Willenbring & von Blanckenburg, 2010). In the study area, no 10 Be flux data are reported by the first two methods, and thus, we adopted the third one, which has previously been applied to large catchments (Rahaman et al., 2017;Wittmann et al., 2015) and very small catchments (<1 km 2 ) (Dannhaus et al., 2018).
We used the output of the fifth-generation GCM, coupled to the aerosol model HAM, which was run separately for 10 Be deposition in the modern (industrial) conditions  and in the early Holocene (preindustrial) conditions . Results of both models were then averaged to represent the mean 10 Be depositional flux in the Holocene (Heikkilä & von Blanckenburg, 2015). The difference between both models is taken as the uncertainty of F 10 Be met . In the Zhuoshui Basin, the average GCMderived F 10 Be met is 1.08 × 10 6 atoms/cm 2 /year, and the difference between both models is small (2.5 × 10 4 atoms/cm 2 /year or 2.3% of F 10 Be met ), indicating a minor climatic impact on 10 Be depositional flux during the Holocene. For obtaining J 10 Be riv , the flux of sediment-carrying reactive 10 Be needs to be known, and we used D insitu as an independent estimate of E in equation (2) due to their similar integration timescales (Willenbring & von Blanckenburg, 2010). In situ 10 Be data (Derrieux et al., 2014) exist near five of our main stem sampling locations (see Figure S6 for locations); hence, we calculated the ratio of exported sedimentary to depositional flux (J 10 Be riv =J 10 Be atm ) for these locations using equations (1) and (2) (Figure 3). Note that for ZS13-6, two D insitu data with highly different values (1.0 and 3.0 mm/year) are nearby; both are used for calculation (see Figure 3), and the resulting difference will be explored later. The derived J   (Table 1) of samples in Figure 3, the estimated F 10 Be met shows a narrow range within the basin, but values are even higher than the GCM-derived one, being 2.42-2.64 × 10 6 atoms/cm 2 /year. The narrow range indicates that the effect of precipitation amount on 10 Be depositional flux would be minor in the studied basin even if it existed. Using the precipitation-based F 10 Be met , the obtained J 10 Be riv =J 10 Be atm ratio would range from 0.09 to 0.50, hence showing a larger flux deficit. These lower ratios may be caused by precipitation records that are not representative of the long-term flux, or by an inaccurate parametrization of the delivery of meteoric 10 Be in this simplified precipitation-based fitting equation. We thus prefer using the GCM-derived F 10 Be met for consistency with our previous studies (Dannhaus et al., 2018;Rahaman et al., 2017;Wittmann et al., 2015). In the future, pending the availability of site-specific F 10 Be met such as from dated soil profiles, more precise D met may be obtained. 5.1.2. 10 Be (meteoric)/ 9 Be-Derived Denudation Rates and Comparison to Other Methods D met are calculated using equation (3) and using a GCM-derivedF 10 Be met of 1.08 × 10 6 atoms/cm 2 /year. Resulting D met generally decrease downstream (Figure 4), ranging from 8.2 mm/year (ZS13-2) to 0.8 mm/year (JRD-S1). Most upstream samples (n = 3) are characterized by D met on the order of 6-8 mm/year. The only exception is ZS17-1 collected closest to the headwaters and in the slate-dominated region, with D met of 3.8 mm/year. Weathered soil profiles in slate lithologies are observed in parts of this catchment where high elevation and high precipitation rates (e.g. >1,800 m and >3.0 m/year;  prevail. Hence, we speculate that the higher [ 10 Be] reac of ZS17-1 (~2-fold of that in the other upstream samples) may partly stem from such local sources that have experienced a longer period of weathering and hence accumulated more 10 Be. Regarding tributaries, both the Chenyoulan and Qingshui Rivers show relatively low D met at 1.0 mm/year. In the downstream floodplain, the denudation rates are mostly uniform at 1.1 ± 0.2 mm/year (n = 4); only the sample ZS17-7 near the river mouth shows a much higher D met of 2.6 mm/year. Considering that the D met of ZS17-7 is quite different from all the other data obtained from a wide time range in the floodplain reaches (bedload collected in 2013 and 2017 and core drilled in 2010), this sample may only record a transient denudation signal and may not be representative of the whole basin. In general, the downstream decreasing trend of D met is explainable by sediment mixing between fast-eroding slate-dominated headwaters and slowly eroding tributaries, which dilutes the high erosion signal from the headwaters in the middle and lower reaches of this basin.
When comparing denudation rates obtained from different methods that integrate over diverse timescales, a key piece of information is the topographic state. As mentioned in section 2.1, the Zhuoshui Basin located in the central Taiwan is likely characterized by topographic steady state. In addition, considering that there is no dependence of denudation rates on time intervals in orogens dominated by fluvial processes (e.g., the Taiwan orogen; Ganti et al., 2016;Sadler & Jerolmack, 2015), denudation proxies intergrating over differnet timescales should be comparable at steady state.
The range of D met in the Zhuoshui River is on the order of 1-10 mm/year, which is consistent with the variation of denudation rates in Taiwan (i.e., several mm/year) over decadal (Dadson et al., 2003), centennial to millennial (Derrieux et al., 2014), and million-year (Hsu et al., 2016) timescales. Especially, the high D met of 4-8 mm/year measured in the upstream slate-dominated Backbone Range is within the same range of exhumation rates (4-8 mm/year) in the core of this mountain belt since~0.5 Ma determined by thermochronology (Hsu et al., 2016;Lee et al., 2006). The agreement between denudation rates from both methods thus supports the use of the GCM-derived F 10 Be met instead of using precipitation-based F 10 Be met . The latter would cause a >2-fold increase of D met for the whole basin and particularly raise the D met in the upper reaches to 10-20 mm/year, much higher than the above-reported exhumation rate.
When comparing to D insitu that integrates over a similar (centennial-millennial) timescale , the sample furthest upstream shows a higher D insitu than those in the lower reaches, showing a trend broadly consistent with the trend of D met (Figure 4a). Nevertheless, two D insitu data within a short distance (upstream distances: 95-107 km) clearly differ (1.0 and 3.0 mm/year) and are generally lower than D met nearby (4.0 ± 0.9 mm/year, ZS13-6). The difference in rates may be attributed to different geomorphic processes recorded in both methods (Portenga et al., 2019) or stochasticity of sediment transport in mountainous rivers (Lupker et al., 2012). At two stations in the middle and lower reaches (upstream distance <80 km), both D insitu and D met agree on the same order of 1-2 mm/year (Figure 4a).
On a shorter (decadal) timescale, modern erosion rates can be derived from suspended load measured by gauging (E gau ; Dadson et al., 2003;Kao & Milliman, 2008). E gau in this study also includes a conservative estimate of bedload transport, that is, 10% of the total sediment load, that is most likely a lower limit for Taiwan rivers (Li et al., 2005;Turowski et al., 2010). The resulting E gau with 4-12 mm/year are distinctly higher than denudation rates derived from cosmogenic nuclides (Figure 4a) for the same river reaches. A reason for this disparity is possibly related to the increasing land use (e.g., road construction; Kao & Liu, 2002) or a higher frequency of typhoon events (Montgomery et al., 2014) in Taiwan during recent decades. If true, E gau may only record a transient signal and could not represent erosion over centennial-millennial timescales probed by cosmogenic nuclides. The higher modern erosion rate is consistent with that in the rivers from tropical highlands affected by agricultural land use (Hewawasam et al., 2003;Vanacker et al., 2007). In contrast, many other small mountainous rivers show higher long-term rates derived from cosmogenic nuclides, explained by the long-term recurrence of large-magnitude sediment-transport events (Covault et al., 2013).

Interplay Between Channel Morphology, Lithology, and Denudation Rates
In orogens with steady-state topography, topography sustains a certain denudation rate to compensate for rock mass accumulation by tectonics (Willett & Brandon, 2002). Hence, a positive correlation between denudation rates and topographic metrics (e.g., channel steepness k sn ) commonly exists (Kirby & Whipple, 2012). However, the coupling between surface denudation and topography can be obscured by spatially variable lithology and climate that can modify the erosional efficiency, and the topography will adjust accordingly based on equation (6). To test the potential control of lithology and climate on denudation processes in the Zhuoshui Basin, we first evaluate if solely a topographic control can explain the variability in denudation rates and then investigate the effects of lithology and climate on the correlation between topography and denudation.

10.1029/2019JF005251
Journal of Geophysical Research: Earth Surface

Inconsistent Spatial Pattern Between Topographic Metrics and Denudation Rates
The correlations between D met and basin-averaged slope, local relief, and k sn are shown in Figure 5. D met is only poorly correlated with these topographic metrics. For example, the correlation coefficient (R 2 ) for basinaveraged slope versus D met is 0.13, and the R 2 for basin-averaged local relief and k sn versus D met are even lower (0.03-0.04, Figure 5). The spatial pattern of basin-averaged k sn along the Zhuoshui main stem (Figure 6a) shows an overall decrease from~187 to~135 km, followed by an increase from~135 to~97 km, with a prominent jump at the confluence with the Danda River in which the local k sn is generally high (Figure 1b), and then followed by a gradual decrease to the downstream floodplain. Hence, the upstream samples exhibit generally low topographic metrics, for example, basin-averaged k sn and local relief (Table 1), but the highest denudation rates (Figure 4). This inconsistent spatial pattern between k sn and denudation rates highlights the complexity of k sn as a denudation indicator.
In order to eliminate the possibility that potential biases in the calculation of k sn contribute to this inconsistency, we consider a range of parameters. One bias arises from the choice of the reference θ value that may not be representative in the studied basin. We thus recalculated basin-averaged k sn using a range of typical θ values (0.35-0.65; Wobus et al., 2006). However, the spatial pattern of k sn remains the same despite of using different θ values (see Figure S1). Another potential bias that can cause abnormal k sn arises from the contribution of flat channels shaped by artificial reservoirs (Figure 4b) or from alluvial regions. After removing the k sn contribution from these flat regions (see Figure S2), the increase of basin-averaged k sn is <5% along the main stem and considered negligible. Hence, we suggest that the basin-averaged k sn , as calculated based on the mentioned assumptions and parameters (see section 3.3) should record a representative pattern of steepness in the Zhuoshui River. Beyond the biases in k sn calculation, the relation between k sn and uplift rates may still be obscured by sediment cover effect and channel width adjustment (Yanites, 2018; Yanites . Both processes commonly reduce the amplitude of k sn increase with increasing uplift rates and hence do not explain the inconsistency between decrease of basin-averaged k sn and increase of D met upstream (e.g., from 90 to 155 km). Overall, our analysis suggests that surface denudation is not solely controlled by topography and lithological and/or climatic factors should play a role.

Impact of Climate and Lithology on the Relation Between Channel Steepness and Denudation Rates
As surface runoff can drive surface erosion (D'Arcy & Whittaker, 2014;Lague et al., 2005), we first evaluate the climatic factor using the correlation between D met and basin-averaged precipitation rate in the Zhuoshui River (Figure 5d). Although the R 2 is relatively high (0.61, Figure 5d), the precipitation rate only decreases by 15% from 3.1 to 2.7 m/year for the Zhuoshui main stem samples (Table 1; Fick & Hijmans, 2017) despite thẽ 8-fold decrease of D met . Considering that denudation rates may be proportional to (precipitation rate) m (m with a common value of 0.4 ± 0.1; D'Arcy & Whittaker, 2014), the minor variation of basin-averaged precipitation rates alone hence cannot explain the variability in D met . Figure 6. Basin-averaged channel steepness (k sn ; a) and areal percentage of each lithological end-member (b) along the main stem. The interval of each dot (a) is~1 km, and one dot includes the topographic information of its whole upstream basin. The low k sn values in the several-km segment of the headwaters may represent a transient feature, caused by basin reorganization and exchange of drainage area on small spatial scales (Chen & Willett, 2016). The areal percentages of the five end-members (b) will always add up to 100%, and the area of Tananao complex and alluvial/terrace deposits is small over basin scale and neglected in (b). The k sn of main stem samples are marked as red open circles in (a) and the sampling locations are marked along the x axis in (b).

Journal of Geophysical Research: Earth Surface
Regarding the effect of lithology on surface denudation in the Taiwan orogen, there is still an ongoing debate. Some studies argued that the difference in lithology is not necessarily related to a change in erodibility or denudation rates (e.g., Fellin et al., 2017;Yanites et al., 2010); other studies suggested that weaker lithologies can increase erosion rates (e.g., Chen et al., 2011;Dadson et al., 2003). Quantifying these effects is complicated, mainly due to two reasons. (1) If the gradients of lithology and rock uplift show a similar trend or the lithology is similar within the studied region, the sole role of lithology is difficult to resolve. (2) Fluvial incision into bedrock occurs by breaking in tension (Sklar & Dietrich, 2001), but rock tensile strength measurements by the Brazilian splitting test (Sklar & Dietrich, 2001) have rarely been carried out in Taiwan (no published data in the study area). On the other hand, rock compressive strength obtained by Schmidt hammer (Selby, 1982) has been measured for each lithological unit in Taiwan Dadson et al., 2003;Siame et al., 2011;Yanites et al., 2010). The compressive strength shows a relatively narrow range for most lithologies of West Taiwan (see the compilation of Figure S3). However, rock compressive strength is not directly linked to fluvial incision and often poorly correlated with tensile strength (Bursztyn et al., 2015). Importantly, the tensile strength of detrital sedimentary rocks can vary by 1 order of magnitude even though the compressive strength is similar (Bursztyn et al., 2015; see Figure S4).
Lithological heterogeneity is large in the Zhuoshui Basin, with rock types varying from schist to sandstone and ages varying from the Late Paleozoic-Mesozoic to Pleistocene (Figure 1b). We illustrate this large heterogeneity in Figure 6b, specifically showing that the areal percentage of both Miocene and Eocene slate in the Backbone Range shows the largest variation compared to other lithological units. We therefore invoke the changes in lithology-related erodibility to explain the inconsistency between channel steepness and denudation rates (Figure 5c), and the workflow is summarized in Figure 7. To quantify the necessary variations in erodibility between lithologies, we applied a revised SPIM that incorporates rock erodibility and river runoff, modified from equation (5) based on Lague (2014) and Sklar and Dietrich (2001): where k * is a constant related to morphological, hydrological, and lithological characteristics, whereas R (runoff, in m/year) and basin-averaged k sn are measurable parameters. Hence, equation (7) needs to be solved for area-weighted e (rock strength index, in MPa) of each lithological unit (see Appendix A for details on model parameterization). k * /e 2 is equivalent to K in equation (5), and thus, the lithology-related erodibility is proportional to e -2 . Note that although the rock strength index here is considered to be mainly controlled by rock tensile strength, the resulting e may include the impact of other lithology-related factors like joint spacing. Smaller spacing may produce more easily removable rock fragments and increase erodibility (Selby, 1982;Whipple et al., 2000).
Determination of acceptable e is based on the following steps (Figure 7). (1) First, we applied a Monte Carlo simulation to randomly generate modeled e of each lithological unit in this basin. Simulations of 10 4 were performed for each unit, and all of them were generated within a reasonable range of tensile strength for  (7)) is described in the box with dashed outline.
(3) We then fitted the D met data set from the Zhuoshui River with the corresponding k Lrsn (10 4 combinations in total); the generated combinations of e of every lithological unit are considered acceptable only when the resulting R 2 value of the correlation is higher than a prescribed threshold (e.g., 0.5). (4) By changing the R 2 threshold within 0.5-0.9 and n within a common range of 1-3 (Harel et al., 2016;Lague, 2014), we can derive acceptable e of each lithological unit for different fitting conditions (i.e., R 2 and n). To test the hypothesis of lithological control, the model output needs to meet three criteria.
(1) There are acceptable values of e that can reconcile variations of multiple parameters, that is, D met , k sn , and areal percentage of each lithological unit and runoff; in other words, the number of solutions for the model should not be zero. (2) The acceptable values of e as a rock property should be largely independent of the choice of topographic parameter (slope exponent n).
( 3) The acceptable values of e should be explainable using independent geological and topographic evidence (Figure 7). When either of these criteria is not met, the hypothesis of lithological control is rejected.
We include four major lithological end-members into this model, including sedimentary rocks of the Western Foothills, (meta-)sandstone of the Hsuehshan Range, Eocene-Oligocene slate of the Hsuehshan Range and Backbone Range, and Miocene slate of the Backbone Range ( Figure 6b). Tananao schist and marble comprise only a negligible area in the whole Zhuoshui Basin (area % <2%) and were therefore not considered. Eocene-Oligocene slate is separated from Miocene slate in terms of their erodibility. The reason is that many thermochronometers are reset in most of the Hsuehshan Range and eastern (Eocene) Backbone Range due to their higher metamorphic grade but are unreset or only partially reset in the western (Miocene) Backbone Range (lower metamorphic grade; Beyssac et al., 2007;Simoes et al., 2012). This difference reflects different burial depth and postdepositional heating and compression history between both units, which can affect the lithology-related erodibility.
Model results for average acceptable values of e of each lithological end-member under diverse fitting conditions (i.e., R 2 and n) are shown in Figure 8a. Note that the relative discrepancy between e of each lithological end-member is largely independent of the choices of R 2 threshold and n and the modeled e of the Miocene slate is always the lowest compared to those of other lithologies. Also, when adding the in situ 10 Be data set (N = 15; see Figure S6 for locations) from West Taiwan (Derrieux et al., 2014), where the lithology and also topographic state are similar to those in the Zhuoshui Basin, model results still show a similar pattern of acceptable values of e ( Figure 8b). Hence, for all kinds of model parameterization, the Miocene slate of the Backbone Range always shows the lowest modeled e and thus the highest lithology-related erodibility. We attribute this high erodibility to its relatively weak compressive strength ( Figure S3), low metamorphic grade (Simoes et al., 2012), and well-developed slaty cleavage, which can become quite fractured during mountain building processes (He, 1986). In support of this, we note that for channel segments draining both Eocene slate and Miocene slate within several kilometers (see Figure 1b or Figure S8 with markers), the local k sn of the Miocene slate is generally lower than that of the Eocene slate. This also suggests a higher erodibility of the Miocene slate according to equation (6) when assuming constant rock uplift within this short distance. In some locations, the absence of a clear contrast in k sn at the lithological boundary can be explained by a gradual transition of lithological properties and/or disturbance by alternations of sandstone and slate (Table S1). Additionally, at the basin scale, the higher k sn in the Danda Basin compared to the upper main stem (Figure 6a) is overall consistent with the lower areal percentage of highly erodible Miocene slate (Figures 8a and 8b) in this tributary (Danda: 18%; upper main stem: 58%).
As expected, when using the average acceptable values of e of each lithological end-member from Figures 8a and 8b, the R 2 for the correlation between denudation rate and k Lrsn increases considerably for D met (R 2 = 0.91) and D met /D insitu (R 2 = 0.78) data sets (Figures 8c and 8d, respectively) compared to the original R 2 of 0.03 (Figure 5c). Hence, a major part of the variation in denudation rates can indeed be explained by variations in k Lrsn . To sum up, given that the modeled results meet the three criteria mentioned above, we propose that the lithological heterogeneity and corresponding variability in erodibility can alter the coupling between denudation rates (recorded in both meteoric and in situ 10 Be) and fluvial channel morphology (recorded in k sn ).

Conclusions
We applied the 10 Be (meteoric)/ 9 Be proxy to investigate the denudation processes in the Zhuoshui Basin, of which~3/5 of the total area is dominated by slates. We examined the controls of topography, climate, and lithology on the spatial variation of denudation rates. The main findings are as follows: Figure 8. Model output from the revised stream power incision model (equation (7)). (a and b) The modeled results of average acceptable rock strength index (e) for different R 2 thresholds (ranging from 0.5 to 0.9) and slope exponents (from 1 to 3, equation (7)), using Zhuoshui D met data set (a, N = 10) and D met and D insitu data sets (Derrieux et al., 2014) in West Taiwan (b, N = 25). We excluded samples in the downstream floodplain and ZS17-1 (with a much lower D met than other upstream samples) to derive representative model results. The y-axis data are normalized to the average acceptable e of the Miocene slate. No solution exists in (b) with a R 2 threshold of 0.9, indicating that factors other than lithology play a role when the data set covers a larger spatial span and is derived from different methods (D met vs. D insitu ). WF = Western Foothills; HR = Hsuehshan Range; BR = Backbone Range; Eo. = Eocene; Mio. = Miocene. (c and d) Examples of the correlation between denudation rates and lithology-and runoff-adjusted k sn (k Lrsn , R θn k sn n /e 2 ) using the modeled e from (a) and (b) with R 2 > 0.5 and n = 1 (see Figure S7 for details of the used e). Data points are colored by area-weighted e. Blue lines and gray dashed lines indicate the fitting trend and the 95% confidence bounds.
1. The range of 10 Be (meteoric)/ 9 Be-derived denudation rates (D met ) in the Zhuoshui River is on the order of 1-10 mm/year, overall consistent with the denudation rates in Taiwan from other methods (sediment gauging, in situ cosmogenic 10 Be, thermochronology) over different respective timescales (from decades to Myr). 2. Along the Zhuoshui River, D met decreases from 4 to 8 mm/year in the slate headwaters to~1 mm/year in the floodplain reaches. This decreasing trend is consistent with the sediment mixing between fasteroding slate-dominated headwaters and slowly eroding tributaries. However, D met is poorly correlated with topographic metrics, especially channel steepness, suggesting that climatic and/or lithological effects complicate the interplay between both metrics in the Zhuoshui Basin. 3. A minor climatic gradient (precipitation rate) in the Zhuoshui Basin fails to explain the large variability in D met . Hence, a revised SPIM that includes the lithology-dependent erodibility is applied. The variations of lithology-and runoff-adjusted k sn can explain the spatial variability in denudation rates as the erodibility of the Miocene slate is predicted to be the highest compared to other lithological units. As such, lithological heterogeneity can indeed alter the coupling between surface denudation and channel morphology.
We conclude that the reasonable denudation rate estimates from the 10 Be (meteoric)/ 9 Be ratio in slate demonstrate the promising prospect of this method to quantify denudation processes in fast-eroding orogens with quartz-poor lithologies. In particular, many small mountainous rivers with high sediment yields are located in high-standing islands of Southeast Asia and Oceania, including global rivers with sediment yields exceeding 10 4 t/km 2 /year (Milliman & Farnsworth, 2011). Considering that the lithology and climate in some of them are similar to those in the Taiwan orogen, the application of the 10 Be (meteoric)/ 9 Be method to quantify their centennial-to millennial-scale denudation processes will be increasingly relevant in the near future.

Appendix A: Parameterization of the Modified Stream Power Incision Model
The constraints on several parameters in equation (7) are explained below. The area-weighted rock strength index e is calculated as where e i is the rock strength index of each lithological end-member and A i the outcrop area of each lithological end-member in a sampling upstream basin. A i of four lithological end-members ( Figure 6b) is derived from the digitalized geological map (Figure 1b). The unknown e i is generated by Monte Carlo simulation (10 4 times). The exponent of e as 2 is based on laboratory abrasion experiments by comparing erosion rate and rock tensile strength (Sklar & Dietrich, 2001), and this factor is considered to provide a reasonable range of erodibilities (Bursztyn et al., 2015).
Regarding runoff (R), an empirical relationship between basin-averaged precipitation rate (Fick & Hijmans, 2017) and measured runoff (Dadson et al., 2003) at gauging stations in Taiwan, is first fitted using a second-order polynomial function. Because gauging records cover inconsistent observation periods between stations and are also insufficient to derive long-term average runoff for some stations, the measured runoff data are relatively scattered for a given precipitation rate ( Figure S5). In this case, the runoff data are binned and averaged at each 0.1-m interval of precipitation rate prior to fitting, to smooth the scatter of gauging records potentially caused by large interannual variability. The final fitting equation, R = −0.025 × P 2 + 0.820 × P (R 2 = 0.63, P represents the precipitation rate; see Figure S5), is then applied to calculate runoff of all sampled upstream basins using precipitation data from Fick and Hijmans (2017). Overall, the impact of runoff is minor in this model given the small climatic gradient in the Zhuoshui Basin (Table 1).
Furthermore, k * as a constant is only used for data fitting and not constrained in the model. The concavity (θ) is fixed at 0.5 as explained in the text, but the slope exponent (n) varies from 1 to 3, considering that the nonlinear k sn -denudation rate correlation can exist in active mountain belts (Kirby & Whipple, 2012).

10.1029/2019JF005251
Journal of Geophysical Research: Earth Surface