Impact of fault damage on eastern Tibet topography

Tectonic deformation can influence spatiotemporal patterns of erosion by changing both base level and the mechanical state of bedrock. Although base-level change and the result-ing erosion are well understood, the impact of tectonic damage on bedrock erodibility has rarely been quantified. Eastern Tibet, a tectonically active region with diverse lithologies and multiple active fault zones, provides a suitable field site to understand how tectonic deformation controls erosion and topography. In this study, we quantified erosion coefficients using the relationship between millennial erosion rates and the corresponding channel steepness. Our work shows a twofold increase in erosion coefficients between basins within 15 km of major faults compared to those beyond 15 km, suggesting that tectonic deformation through seismic shaking and rock damage significantly affects eastern Tibet erosion and topography. This work demonstrates a field-based, quantitative relationship between rock erodibility and fault damage, which has important implications for improving landscape evolution models.


INTRODUCTION
Relationships between topography and millennial erosion rates derived from cosmogenic 10 Be isotopes are typically interpreted to be controlled by tectonics or climate (von Blanckenburg, 2005;Portenga and Bierman, 2011). The topography of the eastern margin of the Tibetan Plateau, which exhibits a >4000 m elevation change over a lateral distance of <100 km ( Fig. 1), has also been attributed to interactions among tectonic uplift, climate, and surface processes (Burchfiel et al., 1995;Ouimet et al., 2009;Kirby and Ouimet, 2011;Scherler et al., 2017). The important role of lithologic compositions and bedrock weathering on erosion and landscape evolution in eastern Tibet was emphasized by Godard et al. (2010) and Gallen et al. (2015). However, impacts of tectonic deformation through seismic shaking and bedrock damage (Faulkner et al., 2010;Huang et al., 2014;Ben-Zion and Zalipian, 2019), which can potentially lead to enhanced erosion (Molnar et al., 2007;Koons et al., 2012;Gallen et al., 2015;Roy et al., 2015Roy et al., , 2016Duvall et al., 2020), have not been quantified in this tectonically active region. Here, we addressed this issue by showing how erosion coefficients, which quantify the relationship between topography and millen-nial-averaged erosion rates, vary as a function of distance to faults in eastern Tibet. To do so, we used cosmogenic 10 Be-derived erosion rates compiled from the literature and our own new measurements. We found a systematic increase in erosion coefficients within 15 km of major faults, and this signal was stronger than that induced by lithologic variation alone.

METHODS
We recalculated previously published 10 Bederived erosion rates (E) from 100 basins across eastern Tibet (Ouimet et al., 2009;Godard et al., 2010;Ansberque et al., 2015) and measured E from an additional 11 basins in the Min Shan area of eastern Tibet. These additional samples came from a tectonically active region underlain mostly by sedimentary rocks ( Fig. 1; see the Supplemental Material 1 and Figure S1 and Tables S1-S2 therein). Please see the Supplemental Material for details about sample preparation and analysis.
We quantified erosion coefficients based on the relationship between E and normalized channel steepness (k sn ). We assumed that (1) river incision into bedrock controls erosion, and (2) detachment-limited bedrock erosion is a function of shear stress or stream power (Howard and Kerby, 1983;Howard et al., 1994;Whipple and Tucker, 1999). Then, we related erosion rate as a function of measurable topographic attributes by: where is the erosion coefficient (a function of climate, channel geometry, and rock erodibility), A is drainage area [L 2 ], S is channel slope [unitless], and m and n are constants. In natural landscapes, local channel slope S can be expressed as, where k s is channel steepness [L 2θ ], and θ is concavity (Flint, 1974;Howard and Kerby, 1983). We extracted channel points with drainage areas larger than 1 km 2 and calculated normalized channel steepness (k sn ) assuming a reference concavity of 0.45, which is consistent with previous studies in this area (Ouimet et al., 2009;Kirby and Ouimet, 2011;Scherler et al., 2017). We calculated basin-averaged k sn by (1) the integral method based on χ (Perron and Royden, 2013;Scherler et al., 2017; see the Supplemental Material) and (2) averaging channel steepness extracted from channel points (Wobus et al., 2006;Ouimet et al., 2009;Kirby and Ouimet, 2011;Scherler et al., 2017). Following Scherler et al. (2017), we used k sn from the integral method to calculate K. By combining Equations 1 and 2, erosion coefficient K can be quantified as ( This relationship holds whether E is in steady state with uplift rate or not. We assumed that sediment production and mixing and grain-size distributions were not significantly affected by local effects. We calculated K assuming a linear relationship between k sn and E (n = 1) following Ouimet (2011) andScherler et al. (2017). We calculated slope as the steepest descent gradient in an 8-cell neighborhood and local relief as the elevation difference between the highest and lowest elevations within a 1-kmradius circular moving window (see the Supplemental Material and Table S3). For sensitivity tests, we also calculated erosion coefficients based on different assumptions of river incision models (e.g., a nonlinear relationship between k sn and E, incorporation of discharge, and different power-law exponents based on stream power or shear stress incision models) or different rate constants calculated from slope and local relief (see the Supplemental Material and Table S4).
We compared K with geologic, climatic, and ecologic factors. We first compiled and quantified the distance to regional-scale (>50 km in length) fault systems (hereafter, major faults) systematically documented by previous studies (e.g., Burchfiel and Chen, 2013). The com-piled faults included both inactive and active (i.e., since the Quaternary) structures, because damage zones of inactive faults could also play a role in shaping present-day topography ( Fig. 1A; see the Supplemental Material). Most of the major faults mapped in this study area are active (86% in total length). Distance to major faults (D mf ) was calculated as the linear horizontal distance from each 90-m-resolution digital elevation model pixel to the nearest fault. We also calculated basin-averaged mean annual precipitation rates (MAP; Bookhagen and Burbank, 2010), normalized-difference-vegetation index (NDVI; Didan, 2015), and peak ground acceleration (PGA) from the 2008 M w 7.9 Wenchuan earthquake (see the Supplemental Material and Figures S2-S3 and Table S3). Basin-averaged values were obtained by averaging all pixel values within studied drainage basins. We also quantified the areal percent of seven lithologic groups for each basin (Table S5) and identified basins dominated by single (defined by >50% of area) lithologies (Table S6).
We performed statistical tests to evaluate controls on K. First, we examined correlations between various controls and erosion coefficients based on different model assumptions and rate constants from topographic metrics (see the Supplemental Material; Table S7). Then, we performed t-tests and F-tests to examine how K was different for subsets of basins divided by D mf , lithologic groups, and basin area (Tables S8-S9). Because large basins may not represent basin-averaged erosion rates due to incomplete sediment mixing (Ouimet et al., 2009;Kirby and Ouimet, 2011), we performed analysis for five basin groups: (1) all basins, (2) basins with A <200 km 2 (hereafter, small basins), and small basins dominated by single lithologies of (3) plutonic, (4) mixed composition sedimentary rocks that included both siliciclastic and carbonate sedimentary rocks, and (5) siliciclastic sedimentary rocks. We focused on three single lithologies that had more than 10 sampled basins. We then used basins that were not significantly affected by the 2008 Wenchuan earthquake for further analysis. A p-value of less than 0.05 was used as the threshold for statistical significance for all tests (see the Supplemental Material).

RESULTS
The 10 Be-derived erosion rates ranged from 0.020 ± 0.002 to 0.76 ± 0.08 mm yr −1 (1σ; Tables S1 and S2). The highest E was from the area near the Huya fault in the Min Shan region (sample ET09: E = 0.76 ± 0.08 mm yr −1 ), and the lowest E was from the southern portion of our study area (sample wbo545: E = 0.020 ± 0.002 mm yr −1 ; Fig. 1A; Ouimet et al., 2009). Basins in the hinterland or low-relief portion of the plateau generally showed lower E values compared to those from the Longmen Shan and Min Shan  (Tables S5 and S6). Averaged D mf ranged from 0.8 to 109.4 km. E showed positive correlations with slope, local relief, and channel steepness, with a wide spread in the data (Fig. 2). Both linear and nonlinear models similarly explained relationships between E and topographic metrics (Table S4).
Both distance to major faults (D mf ) and active faults showed statistically significant correlations with erosion coefficient K and various rate constants for all basins, small basins, and small basins dominated by siliciclastic sedimentary rock ( Fig. 3; Fig. S4; Table S7). The control of D mf on K was observed over a wide range of E and k sn values, which indicates that this result is not from a systematic bias from either E or k sn ( Fig. 2; Fig. S5). Correlations between these various rate constants and MAP, NDVI, and PGA were statistically insignificant or weaker than those with D mf (Fig. S6; Table S7). We found that mean K from basins within a D mf of ∼15 km were ∼2× higher than those with D mf >15 km for all basins and small basins ( Fig. 3B; Table S8). These results were consistent when using erosion coefficients from (1) different river incision model assumptions, (2) different rate constants from topographic metrics, and (3) basins not significantly affected by the 2008 Wenchuan earthquake (Figs. S7-S8; Table S9).
Depending on D mf , similarities between mean K from small basins with different lithologic groups varied. For D mf >15 km, we found no significant differences in mean K for small basins from different lithologies, but there was a statistically significant difference between mean K from siliciclastic sedimentary rocks and mean K from plutonic rocks for D mf within ∼15 km. Lastly, the modeled relation between E and k sn was significantly improved when we considered K varying with D mf , but not with different lithologic groups (F-test, p < 0.05; Table S9).

DISCUSSION AND CONCLUSIONS
Our results indicate that a wide zone (∼15 km) of fault damage significantly impacts erosion coefficients in eastern Tibet, which is consistent with theoretical and numerical studies that emphasize the control of fault damage on erosion and topography (Koons et al., 2012;Roy et al., 2015Roy et al., , 2016. For example, strain generated by faulting may be expressed as newly formed cracks with lengths and density increasing toward faults (Molnar et al., 2007). Bedrock near active or inactive faults that accumulated damage over time is readily detached and transported by erosional processes and significantly influences landscape forms and patterns (Duvall et al., 2020). Considering that most major faults mapped in this area are active, our results are likely derived from impacts of active faults, with potential contributions from inactive faults.
However, our observed ∼15 km width around faults is larger than the typical length scale of fault-damage zones characterized by seismic velocity reduction (e.g., 100-1500 m; Huang et al., 2014) or slow-moving landslide distributions (e.g., 2000 m; Scheingross et al., 2013). This indicates that the impact of tectonic deformation may go beyond localized faultzone cores, which is likely due to combined effects from secondary fault systems and/or coseismic ground shaking (Figs. S2, S3, and S9;Molnar et al., 2007;Faulkner et al., 2010;Pelties et al., 2015;Ben-Zion and Zalipian, 2019). Seismic shaking from large-magnitude earthquakes (>M w 6.5) can enhance erosion by generating widespread coseismic landslides . In fact, most coseismic landslides from the Wenchuan earthquake occurred within ∼30 km of the Beichuan-Yingxiu fault surface rupture (Xu et al., 2014). Based on analysis using these coseismic landslides, Gallen et al. (2015) showed that fractured or weathered bedrock due to high active fault density and precipitation rates may lower rock cohesive strengths (∼65%) near the fault. This is similar to the 71% reduction in overall tensile strength in bedrock inferred from a factor of two increase in erodibility (Sklar and Dietrich, 2001). However, we observed high K values at locations beyond the frontal Longmen Shan, such as the Min Shan, which has complex fault systems consisting of both thrust and strike-slip faults, along which several M w 6 earthquakes have occurred (Xu et al., 2017; see also Fig. 1B; Fig. S3D). This implies that in addition to seismic shaking from earthquakes, long-term plastic deformation of bedrock can be pervasive near major fault zones due to complex fault geometry, kinematics, and tectonic histories.
Our results show that K values from siliciclastic sedimentary rocks are statistically higher than those from plutonic rocks near major faults, which is similar to results obtained via Schmidt hammer rebound measurements (see the Supplemental Material; Figs. S10-S11; Table S10). This may imply that impacts of fault damage on K vary with intrinsic strengths of lithologies. For example, clastic sedimentary rocks often have layered structures and low intrinsic rock strengths (Hoek and Brown, 1997), and so they will be more easily fractured due to static and dynamic stress perturbations near fault zones and will likely produce fine-grained sediments, which are easier to erode or transport (Stock and Montgomery, 1999;Stock et al., 2005). This may explain the pronounced impact of tectonic damage on K in siliciclastic sedimentary rocks from the northeastern region. However, the limited number (n = 12) and coverage (D mf <32 km) of basins dominated by plutonic rocks precludes strong conclusions.
We acknowledge that our fault compilation and data analysis have limitations. First, our compilation may not be complete, because there are likely unrecognized faults due to the long recurrence interval of earthquakes or inaccessibility of field sites. Second, we did not examine impacts from different fault geometries or kinematics due to limited data and complex reactivation histories in this region (e.g., Burchfiel et al., 1995). Certain fault geometries or types (e.g., hanging walls vs. footwalls of thrusts, restraining vs. releasing bends of strike-slip faults) may result in different extents and degrees of fault damage. Third, the quality of fault location and geometry data varied due to the uneven data quality of cited studies. Despite these limitations, our fault compilation provided information on the proximity of faults and allowed us to recognize the impact of fault damage on erosion and topography in eastern Tibet. Our findings imply that the influence of fault damage needs to be considered to adequately model surface processes and examine feedbacks between tectonic deformation and landscape evolution in tectonically active areas (Koons et al., 2012;Roy et al., 2015Roy et al., , 2016. Future studies with more extensive field measurements, detailed information on faults, and increased spatial coverage and number of basins with uniform rock types will help to further clarify the interactions among fault damage, surface processes, and landscape evolution.