# Badland distribution as a marker of rapid tectonic activity

Ci-Jian Yang, Jens M Turowski, Boris Gailleton, Chun-Wei Tseng, Chuan-Min

Chao, Ray Y. Chuang

► To cite this version:

Ci-Jian Yang, Jens M Turowski, Boris Gailleton, Chun-Wei Tseng, Chuan-Min Chao, et al. Bad-land distribution as a marker of rapid tectonic activity. Scientific Reports, 2025, 15 (1), pp.15912. 10.1038/s41598-025-00795-9 . hal-05108345

HAL Id: hal-05108345

https://hal.science/hal-05108345v1

Submitted on 11 Jun 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Badland distribution as a marker of rapid tectonic activity

Ci-Jian Yang<fnref n="1" />✉, Jens M. Turowski<fnref n="2" />, Boris Gailleton³, Chun-Wei Tseng⁴, Chuan-Min Chao<fnref n="1" /> & Ray Y. Chuang<fnref n="1" />

Tectonically driven processes, such as uplift and horizontal advection of rock mass, are important controls on landscape form and dynamics. Available fluvial topographic indicators are intended to be applied to slowly evolving systems over millennial time scales, which may prevent extracting short-term rapid tectonic effects (~1ky). Fast-eroding badlands provide an opportunity to study landscape response to uplift over short time scales. Here, using 1 m airborne LiDAR topographic data in the tectonically active area of southwestern Taiwan, we find that badlands occur within the most rapidly uplifting areas of the landscape, where uplift rates exceed 7 mm y⁻<fnref n="1" />. Furthermore, badlands are predominantly observed within the altitude range that coincides with rapid uplift at 1.3–1.5 kyBP, implying that they formed in response to meso-scale (<100 km<fnref n="2" />) rapid incision of tributaries. We demonstrate how the existing landscape–archive framework, connecting dimensionless erosion rate and relief for hillslopes, can be used on highly erodible badland morphology. As such, badlands can provide archives of uplift history for transient landscapes, which allows to target shorter timescales than previously possible.

Tectonic activity alters landscape form and dynamics and thus mountain building and erosion<fnref n="1" />⁻³. Moreover, it can lead to various natural disasters, thereby increasing the risk to the lives of those residing in affected areas⁴. Constraining spatio-temporal patterns of tectonic uplift help to reveal the mechanisms of regional landscape evolution and supports hazard mapping. Natural landscapes achieve a steady state by adjusting to tectonic uplift over a range of timescales. As such, topography records information on tectonic uplift. For example, fluvial geomorphic indices derived from natural landscapes allow us to obtain tectonic uplift rates and patterns over timescales of hundreds of thousands of years⁵⁻⁹. Yet, the most popular of these indices, channel steepness, does not apply to hillslopes and debris flow channels. Alternatively, the landscape archive framework exploiting hillslope morphology, has been demonstrated to reflect long-term landscape evolution¹⁰⁻¹². With a shorter response time than the river system, this framework may be well suited to constrain the effect of tectonics on topography in low-relief areas with weak lithologies and high erosion rates.

Badlands are fast-evolving landscapes and are especially sensitive to environmental change¹³⁻¹⁵. Badlands typically form due to soft lithologies¹⁶ and climate oscillations¹⁷. Tectonic uplift changes the base level on a regional scale, resulting in incision that propagates upstream through river networks, thereby initiating the expansion of gully and badland systems¹⁸. At the same time, river terraces formed at the onset of the incision can provide temporal reference markers. Thus, once their origin and evolution are understood, fast-eroding landscapes such as badlands, can serve as suitable archives for reconstructing regional uplift histories over short timescales. Here, we suggest that the spatial distribution of badlands reflects the pattern of tectonic activity, and we demonstrate the applicability of the landscape–archive framework to a fast-evolving system.

To elucidate the links between regional tectonic uplift and the distribution of badlands, we study a site located along the active Longchuan fault, on the frontal zone of the Taiwan orogenic belt (Fig. 1a). The Plio-Pleistocene Gutingkeng formation is the primary geological formation. Extending over several kilometers in depth, the main lithology consists of poorly consolidated mudstones, composed of clay (27.6%) and silt (70.6%)¹⁹. Differential erosion and fault activity have been identified as the mechanisms behind the formation of the Longchuan and Qishan faults²⁰, and the associated two regional highlands. These are distinct from the badland landscape evolution caused by short-term tectonic uplift, and we will not discuss them in this study. The Longchuan fault is divided into two segments. The northern segment is a reverse fault in NNE-trending, while the southern segment of the Longchuan fault, also known as the Chegualin fault is characterized as a left-lateral fault with a reverse component²⁰. The timing of tectonic uplift can be constrained since the middle Holocene when multiple river terraces developed in the Erren River, some of which are dated²¹,²². Channel

<footnote n="1">¹Department of Geography, National Taiwan University, No. 1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan (R.O.C.).</footnote>

<footnote n="2">²German Research Centre for Geosciences (GFZ), Telegrafenberg, 14473 Potsdam, Germany. ³Geosciences Rennes,</footnote>

<footnote>University of Rennes 1, 35042 Rennes, France. ⁴Taiwan Forestry Research Institute, No. 53, Nanhai Rd., Zhongzheng</footnote>

<footnote>Dist, Taipei City 10066, Taiwan (R.O.C.). ✉email: cijianyang@ntu.edu.tw</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Study site and background information. (a) Color-shaded relief map showing the distribution of major faults. (b) InSAR results based on ALOS-1 data, showing the LOS displacement, which, shows spatial correlation with mapped badland area. (c) The normalized steepness index of channels ($k_{sn}$, $\theta=0.45$) on a shaded relief map. (d) Grid-averaged annual rainfall map (2000–2010) (source: Central Weather Bureau, https://asrad.pccu.edu.tw/dbar/, last access: 5 May 2025). (e) Regional lithological map (source: Geological map of Taiwan scale 1:50,000, published by Geological Survey and Mining Management Agency, 2013, originally known as the Central Geological Survey, https://www.gsmma.gov.tw/nss/p/index, last access: 5 May 2025).</figcaption></figure>

incision rates into bedrock calculated from the terraces at Yuehshihchieh (Fig. 1a) were 7.6 mm y⁻¹ during ca. 5.7–2.5 ky before present (BP), 5 cm y⁻¹ during 1.5–1.3 kyBP, 2.6 cm y⁻¹ during 1.3–1.0 kyBP and 1.2 cm y⁻¹ since 1 kyBP²¹. The present surface uplift rate (Fig. 1b) has been measured at a maximum of 40 mm y⁻¹ using persistent scatterers interferometry with advanced land observing satellite synthetic aperture radar images from 2007 to 2010 (see details in methods). Neotectonics²¹,²², the development of a diapiric anticline²³,¹ and duplex structure combined with the high-pressure layers within the mudstone formation²⁴ are thought to contribute to this uplift, but the detailed mechanism is still unknown (see details in supplementary). Current topography has been linked to long-term uplift beyond decadal time scales²¹⁻²³. In addition, badlands in southwestern Taiwan have experienced a denudation rate of ~11 mm y⁻¹ at the millennial scale²⁵ and hillslope erosion rates can reach 90–300 mm y⁻¹¹⁴,²⁶.

To connect hillslope morphology with tectonic activity, we used TopoToolbox²⁷ to extract the necessary parameters for dimensional variables²⁸ from 1 m topographic data acquired by airborne LiDAR (see details in methods). To convert altitude into absolute time, we assume that tectonic uplift is uniformly distributed across the spatial domain. In addition, we assume that the stream network rapidly responds to tectonic uplift due to lithological weakness, with a response timescale of decades. Then, all points in the landscape at the same elevation as a particular terrace have the same age as this terrace. Subsequently, we can use hillslope morphology to reconstruct the uplift rate along altitudes of badlands. We use a 1-D landscape evolution model for hillslope morphology to calculate the time-series erosion rate and to invert for uplift rate (see details in methods).

# Results and discussion

## Badlands are concentrated in certain vertical and horizontal areas

Based on a total of 133 million pixels, the average slope gradient of the study site is 22.57 ± 14.37°, and the badland area of 11.9 km² accounts for 9% of the total area. In many landscapes, river bed slope (S) scales with drainage (A) according to a power law, following an equation of the form:

$$S = k_{sn}A^{-\theta} \qquad (1)$$

Here, $\theta$ is a dimensionless exponent known as the concavity index, while $k_{sn}$ is known as the steepness index, which has been interpreted as a proxy for rock uplift rate in steady-state landscapes⁶,⁷. Although the highest decadal uplift rate is ca. 40 mm yr⁻¹ along the hanging wall of the fault, the steepness index $k_{sn}$ is generally low and does not correlate with the uplift pattern (Fig. 1b) at the study site (R²=0.11, $p<0.05$) (Fig. 1c). This indicates that the channels have not reached a steady state slope, possibly due to the influence of ongoing rapid incision into weak bedrock.
Instead, the ratio of badland area positively correlates with surface velocity for surface velocity above 7 mm yr⁻¹ and reaches a maximum value of ~0.68 at a surface velocity of 39 mm yr⁻¹ (Fig. 2a). For median slope gradients above 27°, the ratio of badland area is also positively correlated with the median gradient, and reaches a maximum value of 0.41 at gradients of around 54° (Fig. 2b), and decreases for higher gradients. The ratio of badland area increases from 24 m of elevation to a peak ratio of 0.13 at 40 m and decreases to 0.09 for elevations exceeding 60 m (Fig. 2c).

Our dataset contains a total of 6316 hillslopes, with a mean slope gradient of 0.54 ± 0.22 in tangent value (28.43 ± 12.34°), and drainage areas ranging from 2000 to 33,245 m². The median slope increases from 0.29 to 0.52 at elevations between 23 m and 48 m, it is approximately constant at 0.48 ± 0.01 at intermediate elevations between 55 m and 72 m and decreases to 0.32 again for higher elevations. The median hilltop curvature increases from 0.08 to 0.18 at 41 m and then decreases gradually. Median hillslope length has its highest value of 93 m at an elevation of 25 m and is approximately constant at 64.72 ± 4.69 m for higher values of elevation. Four sets of river terraces can be found in the area. These are at elevations of 25, 35, 50 and 80 m corresponding to ages of 1,

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Correlation of hillslope gradient, LOS displacement, and the ratio of badlands area. (a) Correlation of LOS displacement and ratio of the badland area. (b) Correlation of gradient and ratio of the badland area. (c) Correlation of elevation (0–100 m) and the ratio of the badland area, KT₁ denotes the position of a river terrace. The correlation coefficient is calculated at a significance level of 95%.</figcaption></figure>
1.3, 2, and 6 ky, respectively²¹ (Fig. S4). The largest terrace in the region, named KT₁, was formed as part of the third group and has been associated with ages 2kyBP.

The ratio of badland area positively correlates with surface velocity for surface velocities above 7 mm yr⁻¹ and reaches a maximum value of ~0.68 at a surface velocity of 39 mm yr⁻¹ (Fig. 2a). For median slope gradients above 27°, the ratio of badland area is also positively correlated with the median gradient, and reaches a maximum value of 0.41 at gradients of around 54° (Fig. 2b), and decreases for higher gradients. The ratio of badland area increases from 24 m of elevation to a peak ratio of 0.13 at 40 m and decreases to 0.09 for elevations exceeding 60 m (Fig. 2c).

## Tectonic uplift controls the distribution of Badlands

The formation of badlands has been attributed to a high number of interacting variables, for example, the clay content of the regolith, earthquake-triggered landslides, deforestation, river migration, and land management practices²⁹⁻³². We found that LOS surface velocity derived from SAR estimations positively correlates with the ratio of badland area (Fig. 2a), suggesting uplift related badland evolution. In addition, we found that badlands concentrate at elevations of ca. 40–60 m, corresponding to an age of 1.5–3 kyBP. This may imply that badlands developed during the same rapid uplift period that led to the formation of the KT₁ terrace, the largest river terrace at the study site (cf. Figure 2c). Our interpretation can be summarized in a three-phase conceptual model. The first phase corresponds to the onset of uplift and its initial increase. The main stream incises and forms the first set of terraces. Tributaries respond by quickly incising, forming steep slopes along their channels (Fig. 3a). The second phase corresponds to the rapid increase of uplift until the maximum rate is reached. A second set of terraces forms around the time when the uplift rate is at its maximum. A second tier of steep slopes forms at an elevation below the first tier (Fig. 3b). Both tiers of steep slopes subsequently have developed into badlands by surface erosion and gully incision driven by heavy rainfall. In the third phase, the uplift rate decreases. Still, a third tier of steep slopes forms due to ongoing tributary incision (Fig. 3c). In contrast to previous phases, eroded hillslope material cannot be fully evacuated. It deposits in the upstream reaches of tributaries and turns

<figure><img src="image_5.png" /></figure>

**Fig. 3.** Schematic diagram of the development of mudstone badlands in response to uplift. The dashed line indicates the elevation of the badlands formed with the river terrace; the chronological sequence is distinguished by color. (a) Incision of the main stream (blue arrows) causes terrace formation. Lithologic weakness allows rapid upstream propagation of the incision signal. A first tier of steep slopes forms. (b) Rapidly increasing uplift causes tributary incision, which creates a second tier of steep barren slopes in the headwater regions (pink arrows). (c) Hillslope erosion prevails over decreasing uplift, thus forming sediment-infilled valleys (flat bottomed valley) on the toe of badlands (brown lines). (d) Color-shadow map of sediment-infilled valleys with 1 m spatial resolution.
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Inferential evolution of hillslope morphology. (a) Bedrock incision rate over the 6 kyBP, the elevation and the inferred time are matched based on C terrace dating<sup>21</sup>. (b) Sketch of parameterizations of badland hillslopes. (c) Variation in observed dimensionless erosion rate, *E*\* (Eq. 3, obtained from elevation with respect to dated terraces) with elevation, with colours indicating the binned medians of dimensionless relief *R*\* (Eq. 2). Error bars show the standard error of the mean for binned medians, and solid lines depict binned medians. The simple landscape evolution model (methods) yields modeled *U* (dashed grey line), calculated from Eq. (3), and modeled *E*\* (dashed blue line). (d) Correlation of modeled *E*\* and modeled *U*. Modeled *U* and modeled *E*\* increase between 1 and 1.8 kyBP. Before 1.8 kyBP, uplift decreased, to which erosion showed a delayed response.</figcaption></figure>

them into sediment-infilled valleys (flat bottomed valley) (Fig. 3d). As a consequence of this conceptual model, badlands formed during the same phase have the same age and are located at the same height. Therefore, the topographic features of badlands, such as steep slopes and sharp ridges, should coexist with river terraces at the same elevation. We suggest that these three steps – corresponding to the onset of uplift, its maximum, its decline- and the associated geomorphic evolution provide a general sequence for badland formation in response to rapid tectonic uplift.

## Badlands as transitional landscapes controlled by uplift and erosion

Hillslope erosion rates have been described by a nonlinear diffusion model<sup>11,12</sup>, controlled by hilltop curvature ($C_{HT}$), slope length ($L_H$) and mean hillslope gradient (S). Here, we incorporated the uplift timeline inferred from terrace ages into the nonlinear diffusion model to obtain modeled erosion, $E^*$ and modeled uplift rate, *U* (Fig. 4).

Observed $E^*$ from badland morphology, and modeled $E^*$ increased until ca. 1.5 kyBP (41 m of elevation), when the bedrock incision rate was highest, followed by a decrease (Fig. 4c). Similarly, the modeled *U* driven
by observed E* increases from 1.7 to 6.2 cm y⁻¹ between 1 and 1.9 kyBP, followed by a decrease (Fig. 4d). The modeled U is comparable to the observed uplift rate, supporting the feasibility of reconstructing uplift rates from topography. Within the landscape archive framework¹⁰⁻¹², E* in badland hillslopes can reflect the bedrock incision rate, indicating that badland hillslope erosion is linked to river incision. In line with our conceptual model, we suggest that badland formation is driven by tributary incision in response to uplift, which is reflected in its hillslope form and erosion rate.

In addition, dimensionless relief, R* is influenced by the high slope gradient (Fig. S4) and remains at high values above 0.7 between ca. 2 to 6 kyBP rather than before 2 kyBP (Fig. 4c), suggesting that the slope gradient might be influenced by erosion. Hilltop curvature, C<sub>HT</sub> (Fig. S4) aligns with uplift history (Fig. 4c), indicating that it can effectively record uplift related topography³³. The model predicts that hilltop curvature lags behind relief in its response to changing erosion rates, which makes it possible to distinguish growing from decaying landscapes. Given the correspondence of modeled U and modeled E* (Fig. 4d), we expect that fast-eroding badlands shorten the lag time of topographic response. In addition, they record the rate and timing of topographic response to perturbations during landscape transience.

## Conclusions

Present day surface velocity and hillslope gradient are positively correlated with the ratio of badland area in southwestern Taiwan (Fig. 2a-b), demonstrating the tectonically driven erosion could be expressed in hillslope morphology, mirroring the landscape characteristics in an active orogenic belt³⁴. Elevations with high badland density and modeled erosion rates coincide with the history of rapid uplift at 1.3 to 1.5 kyBP (Fig. 2c), implying that hillslopes have recorded this tectonic event. If badlands can serve as a marker of uplift history, then our observations not only suggest that the existing geomorphic framework applies to erodible badlands but also demonstrate that the integrated change in hillslope morphology with elevation, along with comparable chronological data, can be used to quantify uplift rates between geodetic and geological time scales. Badland landscapes are a microcosm of natural landscapes, observations and numerical simulations based on the rapid hillslope adjustments to tributary incisions caused by tectonic activity can serve as a reference for the landscape evolution of almost all fluvial-connected hillslopes on Earth's surface. Consequently, badlands in this study site serve as an example to address longstanding questions related to the recurrence cycle of geohazards and the evolution of transient landscapes in other rapidly uplifting regions. For example, the concentration of badlands at a particular elevation in the landscape may indicate the occurrence of an earthquake or a change in the activity of a fault.

## Methods

### Measurement of present surface uplift using the PSI method

To characterize present surface velocities (Fig. 1b), we collected 17 advanced land observing satellite synthetic aperture radar (ALOS SAR) images obtained by the Japan aerospace exploration agency (JAXA) from 2007 to 2010, and used the Stanford Method for Persistent Scatterers (StaMPS) software³⁵ to compute line-of-sight (LOS) velocities based on the persistent scatterers interferometry (PSI) methodology. We used fine mode ascending images and excluded the image on March 12, 2010, which contains signals of coseismic and large postseismic displacements of the Jiashian earthquake on March 4, 2010 in southern Taiwan³⁶. The average altitude of the satellite was 691.65 km and the angle between the satellite azimuth and the equator was 98.16°. We used single look complex (SLC) images and chose the image on January 20, 2008 as the master image. We followed the StaMPS procedure³⁵ and selected a coherence value of 0.4 as the threshold to pick PS points. We analyzed a time series of 56 continuously-recorded GPS stations from the GPS lab at the Institute of Earth Science, Academia Sinica³⁷, and calculated their secular velocities. We then projected the 3D GPS velocities onto the LOS direction and compared the projected GPS velocities to the average SAR velocities of all PS points with a 300 m radius around the GPS stations. Finally, we chose the GPS station AKND as a reference point and computed the average SAR velocities concerning the station.

### Extracting parameters from hillslope morphology

To reconstruct the long-term relationship between hillslope morphology and tectonic forcing, we used the landslide inventory database obtained from the Forestry and Nature Conservation Agency (https://www.fore.st.gov.tw/) based on satellite imagery for segment-based interpretation to delineate the extent of barren land. We extracted the slopes for drainage areas above 2000 m², and the mean slope gradient above 5° and manually filtered for artificial structures, i.e., farmland, buildings and roads (Fig. S1). We set the target zone for elevations from 20 m to 100 m to investigate the correspondence between uplift history and hillslope morphology, broadly including all badlands that concentrate at elevations of 50–80 m (Fig. S2 and Table S1). The ratio of badlands was defined as the badlands area divided by the area at the same contour, with values ranging from 0 to 1. For example, we binned the ratio of badlands using bin boundaries at every 1 mm y⁻¹ of LOS displacement and calculated the ratio of badland area within the binned area (Fig. 2a).

Roering et al. (2007)²⁸ suggested a framework of dimensionless variables for analyzing hillslope morphology, consisting of dimensionless relief, R*, dimensionless erosion rate, E*, and dimensionless uplift. The dimensionless relief, R*, is defined by:

$$R^* = \frac{S}{S_c} \qquad (2)$$

Here, S is the slope gradient in tangent value, S<sub>c</sub> = 0.6 is a critical slope value. While this value is below S<sub>c</sub> = 1.2 of Roering et al. (2007)²⁸, the temporal evolution of E* is more crucial for this study. Therefore, in the absence of
better constraints, we choose 0.6 to allow $E^*$ to closely approximate the uplift rate. Dimensionless erosion rate, $E^*$, is defined by:

$$E^* = \frac{2C_{HT}L_H}{S_c} \qquad (3)$$

Here, $C_{HT}$ is the curvature of the hilltop and $L_H$ is the hillslope length.

We binned elevation data using bin boundaries at every 1 m of elevation and calculated the median of $R^*$ and $E^*$ for each bin (Fig. S4). The equations include three topographic features: (i) slope length ($L_H$), is defined as the maximum horizontal distance from the outlet of the slope to the edge of the slope (Figure S3), (ii) slope ($S$), the ratio of a hillslope's height to its length, and (iii) hilltop curvature ($C_{HT}$), the profile curvature of the hillslope divide. All of these values were calculated for each extracted slope by using 1 m topographic data obtained from airborne LiDAR.

## Landscape evolution model for constraining $E^*$ and $U$

We designed a simple 1D numerical model to explore the theoretical response of hillslopes to different uplift forcing and generate time series of $E^*$. Our simulations use the nonlinear diffusion equation originally formulated by Roering (1999)³⁸:

$$\frac{dz}{dt} = \frac{dq_s}{dx} \qquad (4)$$

and

$$q_s = \frac{(K S)}{\left(1 - \left(\frac{S}{S_c}\right)^2\right)} \qquad (5)$$

Where $z$ is the elevation ($L$), $t$ the time ($T$), $qs$ the sediment flux ($L2/T$), $x$ is the distance ($L$), $K$ is the diffusivity coefficient ($L2/T$), $S$ is the local gradient ($L/L$) and $Sc$ the critical slope ($L/L$). This scheme predicts that the long-term evolution of a hillslope tends toward a critical slope $Sc$ with a nonlinear increase of sediment flux as the slope tends toward $Sc$. This equation has been widely used in various forms in numerical simulations³⁹⁴⁰.

The simple nature of our numerical experiments allows us to use a straight-forward centered finite difference scheme to solve the equation explicitly. We use a formulation adapted to run in 1D, as described by Eq. 6 of Perron (2011)³⁹:

$$\frac{(z_x^{t+1} - z_x^t)}{\Delta t} = U + K \left[ \frac{z_x^t}{\left(1 - \left(\frac{z_x^t}{S_c}\right)^2\right)} + 2 \frac{(z_x^{t2} z_{xx}^t)}{S_c^2 \left(1 - \left(\frac{z_x^{t2}}{S_c^2}\right)\right)} \right] \qquad (6)$$

where $U$ is the uplift rate. To set up multiple uplift scenarios, we created an uplift history by linearly interpolating results based on the C dating data. The curve increases from the minimum uplift rate and approaches the maximum uplift rate at 8600 years, then follows the decrease until 10,000 years. After 10,000 years, we set the uplift rate as the minimum uplift rate until the end of the simulation. Note that the maximum uplift rate is set as 10 times the minimum uplift rate, for example, uplift rate = 28 mm/year means the scenario with 28 mm/year of maximum uplift rate and 2.8 mm/year of minimum uplift rate. Therefore, we can keep this time distribution of uplift and scale proportionally to different maximum uplift rates. In total, we conducted 25 simulations with maximum uplift rates ranging from 1 mm/year to 28 mm/year (Fig. S5a). Equation (3) is implemented in Python, using Kwan Lam et al. (2015)⁴¹ JIT compiler to optimize the speed of execution, as the explicit nature of the equation and the rather extreme uplift values tested force small timesteps to keep the scheme stable.

For the initial landscape, we run the model on a trapezoidal hillslope with a low uplift (<0.1 mm y⁻¹) until it reaches a steady state, then we use the low-relief (elevation < 5 m) steady-state hillslope profile with 180 m wide (Fig. S5b) to run our tests obtained by running the model until it reaches stable topography under low-uplift conditions. As such, we can obtain the time series of uplift rates and the corresponding modeled $E^*$ (Fig. S5c). Note that a high diffusivity value of 0.5 was chosen to fit field observation. Our approach remains deliberately simple and intends to explore the theoretical trends in nonlinear hillslopes response to wide ranges of uplifts, as well as time series of modeled $E^*$, rather than reproducing the exact topography observed on the field. Based on the 1D numerical models, the maximum modeled $E^*$'s had a statistically significant correlation with maximum uplift rates (Fig. S5d), which we can use to scale modeled $E^*$ to uplift rate.

To invert the uplift rate from hillslope morphology, we first idealized the observed $E^*$ curve (Fig. 4c) as a Gamma function and determined the shape and scale factors of 1.2 and 98.4, respectively, using a Monte Carlo method. By using the Monte Carlo method, we can flexibly apply the method to any uncertain $E^*$ curve. To keep the shape of the idealized $E^*$ curve and make the maximum $E^* = 1$ for converting to the maximum uplift rate, we normalized the idealized $E^*$ curve with values between 0 and 1. Therefore, we can scale the normalized curve to the modeled uplift rate by the following conversion Eq. (7):

$$modeled\ uplift\ rate = Curve_{normalized} \times (1.392 \times E_{max}^* + 3.397) \qquad (7)$$
## Data availability

Script for the landscape evolution model that supports the findings of this study has been deposited in Figshare https://doi.org/10.6084/m9.figshare.28355798. Other raw datasets are available from the corresponding author on reasonable request.

Received: 6 February 2025; Accepted: 30 April 2025

Published online: 07 May 2025

## References

1. Parker, R. et al. Mass wasting triggered by the 2008 Wenchuan earthquake is greater than orogenic growth. *Nat. Geosci.* **4**, 449–452 (2011).

2. Willett, S. D. Orogeny and orography: the effects of erosion on the structure of mountain belts. *J. Geophys. Res.* **104**, 28957–28981 (1999).

3. Sean, D. et al. Dynamic reorganization of river basins. *Science* **343**, 1248765 (2014).

4. Cattin, R. & Avouac, J. P. Modeling mountain Building and the seismic cycle in the himalaya of Nepal. *J. Geophys. Res.* **105**, 13389–13407 (2000).

5. Howard, A. D. & Kerby, G. Channel changes in Badlands. *Geol. Soc. Am. Bull.* **94**, 739–752 (1983).

6. Kirby, E. & Whipple, K. X. Quantifying differential rock-uplift rates via stream profile analysis. *Geology* **29**, 415–418 (2001).

7. Wobus, C. et al. Tectonics from topography: Procedures, promise, and pitfalls, In Willett, S.D., Hovius, N., Brandon, M.T., and Fisher, D. M. (eds) Tectonics, Climate, and Landscape Evolution, Geological Society of America, Special Paper 398, Penrose Conference Series, 55–74 (2006).

8. DiBiase, R. A., Heimsath, A. M. & Whipple, K. X. Hillslope response to tectonic forcing in threshold landscapes. *Earth Surf. Process. Land.* **37**, 855–865 (2012).

9. Roberts, G. G. & White, N. Estimating uplift rate histories from river profiles using African examples. *J. Geophys. Res.* **115**, (2010).

10. Hurst, M. D., Mudd, S. M., Walcott, R., Attal, M. & Yoo, K. Using hilltop curvature to derive the Spatial distribution of erosion rates. *J. Geophys. Res.* **117**, F02017 (2012).

11. Hurst, M. D., Mudd, S. M., Attal, M. & Hilley, G. Hillslopes record the growth and decay of landscapes. *Science* **341**, 868–871 (2013).

12. Hurst, D., Grieve, S. W. D., Clubb, F. J. & Mudd, S. M. Detection of channel-hillslope coupling along a tectonic gradient. *Earth Planet. Sci. Lett.* **522**, 30–39 (2019).

13. Yang, C. J., Yeh, L. W., Cheng, Y. C., Jen, C. H. & Lin, J. C. Badland erosion and its morphometric features in the tropical monsoon area. *Remote Sens.* **11**, 3051 (2019).

14. Yang, C. J., Turowski, J. M. & Hovius, N. Chang, K.-J. Badland landscape response to individual geomorphic events. *Nat. Commun.* **12**, 4631 (2021).

15. Yang, C. J. et al. Mobile evaporite enhances the cycle of physical–chemical erosion in Badlands. *Earth Surf. Dyn.* **11** (3), 475–486 (2023).

16. Moreno-de las Heras, M. & Gallart F.Lithology controls the regional distribution and morphological diversity of montane mediterranean Badlands in the upper Llobregat basin (Eastern Pyrenees). *Geomorphology* **273**, 107–115 (2016).

17. Calvo-Cases, A. & Harvey, A. M. Morphology and development of selected Badlands in Southeast Spain: implications of climate change. *Earth. Surf. Proc. Land.* **21**, 725–735 (1996).

18. Moreno-de las Heras, M. & Gallart, F. Chap. 2 - The Origin of Badlands, Editor(s): Estela Nadal-Romero, Juan F. Martínez-Murillo, Nikolaus J. Kuhn, Badlands Dynamics in a Context of Global Change, 27–59 (Elsevier, 2018).

19. Lee, D. H., Lin, H. M. & Wu, J. H. The basic properties of mudstone slopes in southwestern Taiwan. *GeoEngineering* **2**(3), 81–95 (2007).

20. Lin, C. W., Yu, C. Y., Hong, G. T. & Chou, P. S. Geological structures of the Tainankaohsiung mudstone area in Southern Taiwan. *Bull. Cent. Geol. Surv. MOEA.* **25**, 143–174 (2012). (in Chinese).

21. Hsieh, M. L. & Knuepfer, P. L. K. Middle-late holocene river terraces in the Erhjen river basin, Southwestern Taiwan-implications of river response to climate change and active tectonic uplift. *Geomorphology* **38**, 337–372 (2001).

22. Hsieh, M. L. & Knuepfer, P. L. K. Synchrony and morphology of holocene river terraces in the Southern Western foothills, Taiwan: A guide to interpreting and correlating erosional river terraces across growing anticlines. *Geol. Soc. Am. Bull.* **55**–74 (2002).

23. Ching, K. E. et al. Rapid deformation rates due to development of diapiric anticline in Southwestern Taiwan from geodetic observations. *Tectonophysics* **692**, 241–251 (2016).

24. Kao, Y. P. Anomalous high deformation rate in mudstone of fold-and-thrust belt in SW Taiwan. Master thesis, Department of Geosciences, National Taiwan University (2017).

25. Hwang, C. E. Suspended sediments of Taiwan rivers and their Geomorphological signifigance. *Bull. Natl. Taiwan. Normal Univ.* **27**, 649–680 (1982).

26. Higuchi, K., Chigira, M. & Lee, D. H. High rates of erosion and rapid weathering in a Plio-Pleistocene mudstone Badland, Taiwan. *Catena* **106**, 68–82 (2013).

27. Schwanghart, W., Scherler, D. & Short Communication TopoToolbox 2-MATLAB-based software for topographic analysis and modeling in Earth surface sciences. *Earth Surf. Dyn.* **2**, 1–7 (2014).

28. Roering, J. J., Perron, J. T. & Kirchner, J. W. Functional relationships between denudation and hillslope form and relief. *Earth Planet. Sci. Lett.* **264** (1–2), 245–258 (2007).

29. Bryan, R. & Yair, A. Perspectives on studies of Badland geomorphology. In *Badland Geomorphology and Piping* (eds Bryan, R. & Yair, A.), 1–13 (Geo Books, 1982).

30. Torri, D., Calzolari, C. & Rodolfi, G. Badlands in changing environments: an introduction. *Catena* **40** (2), 119–125 (2000).

31. Romero Dfaz, A., Sanleandro, P. M., Soriano, A. S., Serrato, F. B. & Faulkner, H. The causes of piping in a set of abandoned agricultural terraces in Southeast Spain. *Catena* **69** (3), 282–293 (2007).

32. Vergari, F., Seta, M. D., Monte, M. D. & Barbieri, M. Badlands denudation hot spots: the role of parental material properties on geomorphic processes in 20 years monitored sites of Southern Tuscany (Italy). *Catena* **106**, 31–41 (2013).

33. Struble, W. T. & Roering, J. J. Hilltop curvature as a proxy for erosion rate: wavelets enable rapid computation and reveal systematic underestimation. *Earth Surf. Dynam.* **9** (5), 1279–1300 (2021).

34. Binnie, S. A., Phillips, W. M., Summerfield, M. A. & Fifield, L. K. Tectonic uplift, threshold hillslopes, and denudation rates in a developing mountain range. *Geology* **35** (8), 743–746 (2007).

35. Hooper, A., Zebker, H., Segall, P. & Kampes A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers. *Geophys. Res. Lett.* **31**(23), (2004).

36. Ching, K. E. et al. Inferred fault geometry and slip distribution of the 2010 Jiashian, Taiwan, earthquake is consistent with a thick-skinned deformation model. *Earth Planet. Sci. Lett.* **31**, 78–86 (2011).

37. GPS lab at the Institute of Earth Science & Sinica, A. (https://gps.earth.sinica.edu.tw) (2022).
38. Roering, J. J., Kirchner, J. W. & Dietrich, W. E. Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology. *Water Resour. Res.* **35** (3), 853–870 (1999).

39. Perron, J. T. Numerical methods for nonlinear hillslope transport laws. *J. Geophys. Res.* **116**, F02021 (2011).

40. Carretier, S., Martinod, P., Reich, M. & Godderis, Y. Modelling sediment clasts transport during landscape evolution. *Earth Surf. Dynam.* **4**, 237–251 (2016).

41. Siu Kwan, L., Antoine, P. & Stanley, S. Numba: a LLVM-based Python JIT compiler. In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC (LLVM '15). Association for Computing Machinery, New York, NY, USA, Article 7, 1–6 (2015).

## Acknowledgements

We express our gratitude to the anonymous reviewers for their constructive suggestions, which have greatly improved this study. We apologize for not being able to list each individual's name. This study was supported by grants from the National Science and Technology Council of Taiwan awarded to C.-J. Yang (MOST 110-2917-I-564-009 and NSTC 112-2116-M-002-033-MY3).

## Author contributions

C.-J.Y. designed the study, data analysis, and modelling. J.M.T. and B. G. designed modelling. C.-W.T. provided the high-resolution spatial data. C.-M. C. and R. Y. C. analyze surface velocities driven from synthetic aperture radar. C.-J.Y wrote the paper with input of all authors.

## Declarations

### Competing interests

The authors declare no competing interests.

### Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-025-00795-9.

Correspondence and requests for materials should be addressed to C.-J.Y.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025