# Deformed alluvial terraces record an excess of slip over the last few centuries on the Himalayan Topographic Frontal Thrust of central Bhutan

Maxime Gautier, Matthieu Ferry, Stéphanie Gautier, Romain Le Roux-Mallouf, Dowchu Drukpa, Rodolphe Cattin

► To cite this version:

Maxime Gautier, Matthieu Ferry, Stéphanie Gautier, Romain Le Roux-Mallouf, Dowchu Drukpa, et al.. Deformed alluvial terraces record an excess of slip over the last few centuries on the Himalayan Topographic Frontal Thrust of central Bhutan. Communications Earth & Environment, 2024, 5, pp.590. 10.1038/s43247-024-01759-z. hal-05459418

HAL Id: hal-05459418

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

Submitted on 15 Jan 2026

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Deformed alluvial terraces record an excess of slip over the last few centuries on the Himalayan Topographic Frontal Thrust of central Bhutan

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

Check for updates

Maxime Gautier<fnref n="1" />,⁴, Matthieu Ferry <fnref n="1" />,⁴, Stéphanie Gautier<fnref n="1" />,⁴, Romain Le Roux-Mallouf²,
Dowchu Drukpa ¹³ & Rodolphe Cattin <fnref n="1" />,⁴ ☑

Deformed alluvial terraces are ubiquitous markers of a fault's recent activity and may help assess its slip rate and associated seismic hazard. They are often considered as a nearly flat surface translated and rotated along a planar or listric fault. The present study challenges these assumptions by revealing uneven terrace treads and verticalization of the Topographic Frontal Thrust (TFT) in south-central Bhutan. We model this finding as combined variability in both the aggradation and geometry of the TFT. We estimate a Holocene slip rate of 19.6 ± 4.1 mm.yr⁻<fnref n="1" />, which confirms that the TFT accommodates most of the shortening across the range. Contrary to previous studies, we find an excess of slip over the last few centuries, which implies a lower seismic hazard. These results highlight the importance of considering the non-planar component in terrace shape, shallow abrupt changes in fault geometry, and aggradation in future morphotectonic studies worldwide.

Nearly ubiquitous along fluvial valleys, river terraces are major geomorphological archives of climate and tectonic processes (e.g. ref. 1). The extensive related scientific literature underscores the importance of these landforms. Numerous studies focus on terrace formation through field observations to survey their geometry and collect samples for determining their age of deposition or abandonment (e.g. ref. 2), and digital elevation model analysis to map these landforms<fnref n="3" />. In parallel, a whole sub-field of geomorphological studies carries out experimental approaches to explore the mechanisms of channel response and resulting terrace geometry to changes in hydrological regimes (e.g. ref. 4), and numerical modeling to study their emplacement through climatic modulation (e.g. ref. 5).

Because they represent surfaces of potentially known age with a relatively simple initial geometry, alluvial terraces affected by active faults are often used to quantify late Quaternary tectonic activity and the associated seismic hazard (e.g. refs. 6–9). In the most straightforward cases, the preserved terraces are approximately planar landforms uplifted by thrusting at depth. To estimate an uplift rate, Quaternary dating methods (generally radiocarbon, luminescence, or cosmogenic exposure) are combined with topographic profiles surveyed over the terrace treads, for which elevation is compared to the current river bed considered at constant equilibrium. This

rate is then converted to fault slip by assuming a rigid block model with a 2D fault geometry defined by a constant dip angle estimated from surface observations⁸,¹⁰⁻¹⁴. However, when a terrace is uplifted and deformed by the activity of a thrust fault, the information is sometimes more difficult to extract as the terrace tread may exhibit a non-planar shape. Converting observed surface terrace tread elevation into fault slip requires more complex modeling accounting for fault dip variations at depth and incremental slip controlling progressive deformation (e.g. refs. 15–19). These kinematic models have successfully been applied to folded and tilted river terraces above listric thrust faults, representing a subset of purely concave structures (see ref. 16 for a comprehensive and detailed example). A more general framework is still needed to thoroughly explore geometrical and mechanical parameters controlling the unevenness of alluvial terrace treads.

In central Bhutan, the tectonic activity is mainly controlled by the Himalayan Topographic Frontal Thrust (TFT), the most recent expression of the thrust sequence that accommodated the India-Eurasia collision over geological timescales. In the past millennium, at least two great seismic ruptures occurred along this fault: an Mw 7.5–8.5 earthquake in 1714 CE and an Mw >8.5 earthquake that produced ca. 8 m of uplift during medieval times²⁰. Alluvial terraces, formed by monsoonal episodes and earthquakes,

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

e-mail: rodolphe.cattin@umontpellier.fr

<footnote n="1">¹Géosciences Montpellier, Université de Montpellier, CNRS, Montpellier, France. ²EDF-DIPNN-DI-TEGG, Pôle aléa sismique, EDF, Aix-en-Provence, France.</footnote>

<footnote n="3">³Department of Geology and Mines, Thimphu, Bhutan. ⁴These authors contributed equally: Maxime Gautier, Matthieu Ferry, Stéphanie Gautier, Rodolphe Cattin.</footnote>
are preserved along the whole Himalayan front and especially well in the Sarpang area of central Bhutan (Fig. 1). The study of Berthet et al.¹² on the main alluvial terrace suggests a Holocene uplift rate of 8.8 ± 2.1 mm.yr⁻¹. These authors consider a rigid block model with a fault dip of 25° ± 5° and obtain a thrusting rate of 22.6 ± 9.2 mm.yr⁻¹. This result presents large overall uncertainties stemming from poorly constrained terrace deformation, fault structure, and fault dip. Consequently, the role played by this frontal fault in accommodating shortening across the Himalayan range remains unresolved, although it is crucial information for estimating seismic hazards along the Bhutan Himalayas.

Here, we propose an upgraded model for interpreting deformed allu-
vial terraces such as these observed in the Himalayas (Fig. 2). This model,
while similar to kinematic models in its ability to replicate the successive

deformations of terrace treads caused by seismic or creep events along a listric thrust fault, surpasses the intrinsic limitations of such models by taking into account abrupt changes in fault dip and curvature, mechanical properties of superficial materials, and aggradation. Our approach relies on a refined geophysical image of the TFT beyond a depth of 5 m and a detailed mapping of the unevenness of terrace treads through field observations and a very high-resolution digital elevation model (see “Methods” section). Our modeling strategy consists of two main stages (Fig. 3). Firstly, we gather a priori information to establish a plausible range of model parameter values. Field observations and geophysical measurements aid in defining the geometry of the thrust fault at surface and depth, as well as the thickness of aggraded sediments in the current riverbed. Secondly, we adopt a stochastic approach to explore the model’s parameter space. Each model’s likelihood is

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

**Fig. 1** | Geomorphological setting of the Sarpang area within the Himalayas. a 3D oblique view of Pleiades-derived Digital Surface Model (DSM) at 1 m resolution (see Method section for details) with a vertical exaggeration of 2.5. The Topographic Frontal Thrust (fault surface trace, black line) separates the undeformed alluvium of the Indian plain from a set of abandoned terrace treads (T1–T7) that exhibit mean elevations above the present stream belt (blue lines draining T0) ranging from ca. 8 m to ca. 64 m. Gray lines are topographic profiles P0-P8 extracted from Pleiades-

derived DSM (x-origin next to profile label). The white triangle represents a soil pit sampled by Berthet et al.¹² to date the deposition of T6. Yellow diamonds are locations of paleoseismic trenches from Le Roux-Mallouf et al.²⁰. The top surface of T6 is colored according to elevations, as shown by the color scale. Green line shows the location of the electrical resistivity profile. The viewpoint of Fig. 5a is also given.

b Transverse west-east topographic profile P8 showing the sequence of fluvial ter-
races of the Sarpang river.
<figure><img src="image_6.png" /><figcaption>(b) Upgraded model</figcaption></figure>

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

**Fig. 2** | Simplified sketch illustrating the various processes impacting the shape of alluvial terrace tread in the Sarpang area. **a** The conventional approach used by Berthet et al.¹² views terraces as markers of tectonic activity associated with movement along planar fault surfaces with constant dips. The alluvial plain deformed by these faults is a flat surface over which a channelized river flow. The resulting terrace tread exhibits a smooth surface, and its geometry can be associated with the constant uplift of a rigid block. **b** Geological and geophysical observations indicate that more than this simple approach may be required. Alluvial deposits, even in active river beds, can vary in height by several meters, with geomorphic markers such as natural levees, point bars, banks, or pools. We propose an upgraded model to interpret the observed unevenness of terrace tread, accounting for a complex fault geometry, as suggested by lateral variations of the scarp toe or imaged by geophysical investigations and considering variations in aggradation.

assessed by comparing the simulated topography with the topography of terrace treads extracted from a very high-resolution digital elevation model. This approach yields a posteriori insight into the thrust fault geometry and initial terrace elevation. The obtained cumulative displacement is then combined with age constraints to derive fault slip velocity. We apply this approach to the alluvial terraces observed in the Sarpang area and derive the slip rate distribution along the TFT (Fig. 1). A comparison of the slip rate distribution we obtained for the different terrace generations with the geodetic shortening velocity measured in central Bhutan leads us to re-evaluate the accommodation of deformation in this region, particularly during the Holocene.

# Result

## Geometry and composition of the Sarpang terraces

To document in detail the geometry of the Sarpang terraces, we analyze the digital elevation model (DEM) and extract topographic profiles (Figs. 1 and 4). We identify a staircase of eight main terraces, T0 to T7, from

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

**Fig. 3** | Workflow describing the two main steps of our approach. The approach includes¹ model setup and a constraint on model parameters through geophysical imagery and field observations, and² stochastic approach for which the simulated topography is compared with the topography of terrace treads extracted from a very high-resolution digital elevation model. This approach yields posterior information on the fault geometry and height of pre-deformed alluvial deposit. The obtained cumulative displacement is combined with age constraints to derive fault slip velocity.

lower to higher with respect to the present stream with relative elevations of 8–8.5 m (T1), 11–11.5 m (T2), 12.5–13 m (T2b), 22.5–23 m (T3), 34.5–35 m (T4), 37.5–38 m (T5 or possibly T4b), 58–60 m (T6), and 62–64 m (T7) (Fig. 1). All terraces are apparently flat, exhibiting a general gentle slope of ca. 1% to the south (direction of flow from the main drainage) and are bounded by hillslopes to the north and by steep tectonic scarps to the south along the TFT (Fig. 1a). Terrace treads show continuous deposits of polygenic well-rounded pebbles and blocks. Terraces T2 and T6 (Figs. 1 and 4) have been extensively studied and their age has been determined by previous authors¹²,²⁰,²¹. They are the focus of the present work. Though well documented, other terraces (T2b, T3, T4, T5 and T7) exhibit unfavorable geometry, extensive reworking/erosion, or densely vegetated scarps and are less suited to accurately record and document deformation.

Furthermore, the active river bed T0 presents significant rugosity with bars. It swells reaching ca. 4 m elevation above the present stream level (as of March 8, 2013, the acquisition date of Pleiades images). In the Sarpang area, the climate is dominated by the monsoon season that starts in April and ends in October, with a peak of approximately 1000 mm of rainfall in July (Climate Data https://en.climate-data.org/asia/bhutan/sarpang-district/sarpang-707053). The rest of the year experiences minimal rainfall. Hence, sediment transport and deposition occur during and after the monsoon peak and are followed by seasonal downcutting during the dry
<figure><img src="image_9.png" /><figcaption>**Fig. 4** | Bulge width and height measurements from topographic profiles. The red line is a linear approximation of the present stream (profile P0) with a constant slope of 0.9°. The northern section of each profile is fitted to that trend. The part of the profile above the trend along the southern section is defined as the bulge. Dashed</figcaption></figure>

lines indicate how the dimensions of the bulge are defined along with values and uncertainties. For profiles P0 to P7, x-axis is the distance from the scarp toe in meters, y-axis is the elevation in meters.

season, which is the signal recorded by our Pleiades-derived DEM. Hence, the approximately 4-m-high seasonal water level variability should be considered and not confused by tectonic-driven downcutting. As expected, the most significant variability is observed along the most active channels and decreases in their vicinity. On a broader scale, profile P0 (Fig. 4) shows that T0 is an overall planar surface with a constant slope of 0.9° in the North-South direction of average flow.

Terrace T2 (Figs. 1 and 4) exhibits an identical slope to T0 but at higher elevations by 5 ± 1 m and forms a wide bulge (370 ± 20 m) whose summit rises to 9.5 ± 1 m relative elevation. The southern tip of that 4.5-meters-high bulge corresponds to a ca. 4-meter-high tectonic scarp (Fig. 1 upper left panel) documented by Berthet et al.¹², where Le Roux-Mallouf et al.²⁰ observe the emergence of the TFT. The latter authors show that T2 comprises a 3-meter-thick deposit of pebbles to cobbles with a sandy matrix overlain with a 1-m-thick sand and silt layer that sits directly over the bedrock (Lesser Himalaya sequence). There, the fault contact is horizontal

as it intersects the topography, then steepens to ca. 12° over a few meters along the dip, then to ca. 27° along the T2/bedrock contact as the fault enters active bed sediments (T0).

Terrace T6 is a well-developed regional surface that extends continuously over ca. 1.5 km². Dense vegetation and steep scarps hinder the observation of natural exposures in the field. However, soil pits excavated by Berthet et al.¹² show that alluvial deposits are at least 3 m thick and composed of pebbles to blocks within a sandy matrix. Profiles P2 to P7 (Fig. 4), extracted along the assumed flow direction, consistently display a flat surface along the northern half of the tread with a southward slope of 0.9° and a bulge along the southern half with a width of ca. 530 and a height of 13.8 ± 2 m above the linear trend of the northern section (Fig. 4). We interpret the linear sections of the profiles to reflect the original geometry of the stream profile (as documented by T0/P0 in Fig. 4) displaced by translation along a linear fault section and the bulge to record folding over a curved fault section. It should be noted that the bulges on profiles P6 and P7
<caption>**Fig. 5** | Multi-scale multi-method characterization of the recent deformation at the Sarpang river site.</caption>

a Field photograph (see location Fig. 1) shows partially exhumed bedrock (green line is top of Lower Siwaliks) overlain with ca. 2.5-m-thick alluvium from T2. T2 is uplifted by ca. 4 m along the emergence of the Topographic Frontal Thrust (black line)¹². Note T2 is incised by ca. 1.5 m for the present stream level. This outcrop was later cleaned and documented in detail for a paleoseimological study by Le Roux-Mallouf et al.²⁰, who measured a near-horizontal fault dip close to the surface steepening to 27° at the base of the outcrop. b General topography of T0, T2, and T6 terraces at the Sarpang Chu site. Profiles are extracted from the processed photo-grammetric DSM (see Methods for details). Field observations (panel a) are reported to scale and correspond to the scarp affecting T2 visible in the profile. The lower insert is the ERT profile (panel c) reported to scale. c Merged multi-scale ERT profile surveyed at the same location. Black inverted triangles indicate the location of electrodes for the three original profiles. Thick black lines indicate the fault's geometry at depth tested by modeling, i.e., α = 20°–30° for the shallow section and β = 50°–80° for the deeper section.

South

North

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

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

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

(Fig. 4) appear reduced in amplitude (ca. 12 m versus 14-15 m for P2 to P5). Furthermore, the planar sections of profiles P2-P5 display elevations that are systematically 6–7 m lower than the planar sections of profiles P6-P7. Since the origin of these differences is unclear (e.g., lateral variation, erosion, incomplete DEM due to very dense vegetation), they will be modeled separately hereafter.

## Near-surface TFT imaging

We complement the geomorphological study with geological information on the TFT to unravel the effect of tectonic uplift and aggradation on these terrace features. Structural measurements in the Sarpang area give a dip of approximately 30° towards the North²², which is consistent with the dips of 20°–30° towards the North observed below T2¹²,²⁰ (Fig. 5a). We also use near-surface information from shallow geophysical surveys conducted along the Sarpang River near the scarp toe of T2 (Fig. 1). We combine two complementary methods¹: a nested electrical resistivity tomography (ERT) providing a resistivity model and its reliability²³,²⁴ and a stochastic approach to account for resistivity contrast on both sides of the TFT²¹.

The ERT image and the analysis of the misfit distributions lead to a robust and accurate characterization of the fault geometry down to a depth of ca. 40 m (Fig. 5c and Supplementary Fig. 1). The set of most probable models derived from the stochastic approach explains the main resistivity pattern of the ERT images. However, the histograms point out that the thickness of the different layers remains poorly resolved. The distribution of resistivity values delimits three distinct bodies. South of the emergence of the TFT, a homogeneous unit at least 10 m in thickness (local maximum penetration depth) displaying high resistivity values of 1000–6000 Ohm.m that reflect the alluvial deposits of the Indian Plain observed at the surface. North of the TFT, a shallow body ca. 5 m in thickness to the north thins down to approximately 3 m near the fault trace displays resistivity values of ca. 1000 Ohm.m and matches the alluvial deposits from T1. Underneath, a

massive unit is imaged down to ca. 40 m (maximum penetration depth), displaying resistivity values ranging from 10 to 100 Ohm.m (Fig. 5c). We associate this signature with the local exhumed bedrock composed of weathered Lesser Himalaya phyllites²² that contain a high proportion of water-saturated clay minerals²⁵ known to be electrically conductive²⁶,²⁷. A sharp and well-localized contrast zone separates this latter unit from the previous one. It suggests that the tectonic contact between recent alluvium and the overly thrusting bedrock is shallow near the fault and steepens further north.

Our results suggest that the fault dip angle observed at the surface is only valid at the very shallow depth. Indeed, at a depth of ca. 5 m, our findings reveal an abrupt and significant steepening of the fault from 20°–30° to 70°–90° (Fig. 5c and Supplementary Fig. 1). This pattern is coherent with the flat and listric-ramp geometry proposed for the TFT from gravity measurements²¹, suggesting that the fault's geometry changes dramatically with depth.

## Modeling of the tread terraces unevenness

We examine the unevenness of T2 and T6 through three sets of topographic profiles from west to east: P1, P2-P5, and P6-P7. Based on elastic dislocation full-space modeling, our approach compares an observed elevation profile to profiles produced by different fault geometries, cumulative displacement, number of seismic events, and height of pre-deformed alluvial deposit, which are estimated a priori from field observations and geophysical investigations (see “Methods” section).

The modeled topographic profiles agree with the elevation of the terrace treads, accounting for the total vertical offset and the bulge's amplitude, wavelength, and location observed for T2 and T6 (Fig. 6). The results gathered for profile P1 across terrace T2 (Fig. 6d) reveal a concave geometry with a dip decreasing with increasing depth. This model does not satisfy the simple flat-and-ramp structure generally hypothesized for
<figure><img src="image_13.png" /></figure>

**Fig. 6** | Comparison between observed and best-fit modeled terrace tread elevation. a Gray circles represent fault-perpendicular topographic profiles P6 and P7 extracted from Pleiades-derived DSM across the deformed terrace T6. Blue dashed lines show the vertical displacement calculated with a rigid block model for a fault displacement ranging from 119 m to 159 m as proposed by Berthet et al.¹². The blue solid line is the best-fit surface deformation obtained with the proposed modeling approach for a cumulative displacement of 138 m and taking into account aggradation and fault dip angle variation (see white star marker location on Supplementary Fig. 2 for the value of model parameters). b Same as a for topographic profiles from P2 to P5. The best-fit surface (green solid line) is obtained with a fault displacement of 136 m (see white star marker location on Fig. 8 for the value of model parameters). c Same as a for profile P1 across the terrace T2. The best-fit surface (red solid line) is obtained with a fault displacement of 12.9 m (see white star marker location on Fig. 7 for the value of model parameters). d Simplified sketch showing the geometry of the two tested models. The dashed line gives the rigid block model geometry that assumes a constant fault dip angle defined by the slope observed at the surface. The solid line represents the proposed model based on an elastic dislocation full-space approach with four north-dipping segments, each with its dip angle and length (see text and Figs. 7 and 9 for details about trade-offs and distribution of these parameters). For simplicity, we represent the fault geometry associated with profile P1, knowing it varies laterally and sediment thicknesses differ between T2 and T6.

the Main Himalayan Thrust but is consistent with a flat-and-listric geometry where the listric section would be rooted into a subhorizontal level at a depth of at least 300 m. They also suggest a cumulative fault slip of 15.8 m and a height of pre-deformed alluvial deposit of ca. 3 m (Fig. 7). Despite the small cumulative displacements resulting in a low signal-to-noise ratio, the well-characterized nature of T2 enables us to test the robustness of our results. As depicted in Fig. 5a, field observations suggest a vertical tectonic offset of ca. 4 m close to the fault, consistent with a 10-20 m slip along a fault with a dip of 20°–30°. The relative location of the fault with respect to the modeled topography also agrees with field observations¹²,²⁰. Lastly, T2 observed in the hanging wall of the

fault features an alluvial deposit of ca. 3 m, aligning with our result (Fig. 5a). In parallel, our inversion of swath profiles between P2–P5 and P6–P7 across T6 gives a significantly larger amount of cumulative slip, reaching up to 138 m (Fig. 6a, b). Over the investigated depth, the modeled TFT geometry is similar, with a steeper and straighter flat-ramp-flat geometry (Fig. 6d). The assessed thickness of pre-existing alluvial deposits ranges between 4 and 10 m (Fig. 8 and Supplementary Fig. 2).

The cumulative slip can also be estimated from Eq. (1) presented in “Methods” section for the flat threads approach. Assuming a constant dip angle of 25°, this rigid block model requires two values for the cumulative slip to account for the tread terrace unevenness. The total slip from west to east for P1, P2–P5, and P6–P7 is then estimated to be 2–13 m, 105–145 m, and 119–159 m, respectively (Fig. 6). This simple approach gives results consistent with our findings, but may miss part of the geomorphological signal, leading to misinterpretation. For instance, with the rigid block model, the lateral variation and the associated 6–7 m difference in elevation between profiles P2–P5 and P6–P7 are interpreted as a 14 m along-strike difference in slip (Fig. 6a, b). In contrast, our approach may explain this apparent eastward increase in cumulative uplift from a combined effect of four model parameters: variations in slip (138 m vs. 136 m), thickness of pre-existing alluvial deposit (7.5 m vs. 3.5 m), fault dip angle (27° vs. 23.5°) and distance between the scarp toe and the fault (110 m vs. 13 m) (Figs. 6, 8 and Supplementary Fig. 2). Our approach highlights the potential role of a lateral variation in aggradation (consistent with the distance to the Sarpang River) as well as along-strike and across-strike changes of the TFT geometry, rather than attributing the along-strike variations observed on the T6 terrace solely to an abrupt eastward increase in slip.

# Discussion

## Model parameters trade-off and distribution

In the previous paragraph, we have only presented the results associated with the best-fitting models. The two million calculated models also allow us to characterize the interdependence among the model parameters and their distribution associated with the most likely models. The inversion of topographic profile P1 reveals substantial trade-offs between the fault dip angle, the total amount of slip, and the thickness of alluvial deposits (Fig. 7). Overestimating (underestimating) the dip angle of the deepest fault segment by a few degrees or underestimating (overestimating) the height of alluvial deposits by a few meters can lead to an overestimation (underestimation) of the fault slip by several meters. Our results suggest a linear relationship between the thickness of pre-existing alluvial deposits *H* and the cumulative slip *U*<sub>cumulative</sub>, with the sum of these two parameters equaling 17–19 m. The analysis of profiles P2–P5 (Fig. 8) and P6–P7 (Supplementary Fig. 2) confirms this finding, with a clear linear relationship between *H* and *U*<sub>cumulative</sub>. However, it shows that the relative effect of aggradation is less significant as the height of the terrace increases. Therefore, the greater the relative elevation of the terraces, the better the constraint on fault geometry and cumulative slip. Since footwall terraces can be poorly preserved, this result underlines the importance of correctly deciphering the morphological signal of the unevenness of terrace treads in the hanging wall when the relative elevation of preserved terraces is comparable to the rugosity of the active river bed.

As the two million scenarios tested are randomly drawn with equiprobability, the study of the most plausible models yields a probability density function (pdf) for each model parameter (Fig. 9). The densities obtained by considering the thousand best-fitting models show a flat distribution for the number of events *N*, suggesting that this parameter has no significant effect on the morphology of terrace treads. This result could be due to the study's failure to consider diffusion processes. Although the slight variations in slope on the terrace tread surface justified neglect of the erosion process, this assumption becomes invalid to model the degradation of terrace risers, which can be used to constrain fault slip history and the associated number of events²⁸.
<figure><img src="image_14.png" /><figcaption>**Fig. 7** | Estimated trade-offs between model parameters for terrace T2. Estimated trade-offs between the geometry of the steepest fault segment (dip angles $γ$ and $δ$), the initial thickness of pre-deformed alluvial deposit *H* and the cumulative displacement $U_{cumulative}$ derived from topographic profile P1 across the terrace T2. This scatter plot is obtained from the 2,000,000 models stochastically explored. Circles are</figcaption></figure>

colored according to model likelihood as shown by the color scale. The likelihood ranges from 0 (poor agreement with the observations) to 1 (perfect fit). White star shows the best-fitting model.

Similarly, the pdf associated with the length $L_1$ and $L_2$ and dip $α$ and $β$ of the most superficial fault segments do not provide any posterior information for these parameters, which are already well-constrained by geological and geophysical observations. The geometry of the deepest segments varies significantly between profiles P1 and P2-P7 (Fig. 9). The results show a listric fault geometry for P1 and a near-constantly dipping fault of 20°–25° and 23°–28° for P2–P5 and P6–P7, respectively. The densities associated with cumulative displacement exhibit a normal distribution with a mean of 16 m, 134 m, and 136 m and a standard deviation of 4.6 m, 9 m, and 8 m for P1, P2–P5, and P6–P7, respectively. The obtained distributions confirm little or no west-east variation in cumulative slip along the TFT under the T6 terrace, which means that the amount of slip is relatively consistent from the west to the east, as previously mentioned in the “Result” section for the best-fitting model.

## Fault slip rate and seismic hazard implications

Calculating the slip rate is critical for assessing fault activity and the associated seismic hazard. Here, we reevaluate the slip rate distributions from the

available age and the obtained cumulative displacement distributions presented above.

The ages of terraces T2 and T6 were originally determined by Berthet et al.¹² (Fig. 10a). They excavated a 70-cm-deep pit into the top deposits of T2 (sand and silt unit) and sampled detrital charcoal lumps that were subsequently radiocarbon-dated to 1550 ± 70 CE (see Fig. 1 for the location of the westernmost trench). These dates coincide with the timing of paleoseismic events affecting the lower T2 units deposited directly over a strath surface on top of the Lesser Himalaya bedrock²⁰. The former also excavated a 3-m-deep pit within the top surface of T6 over a flat section likely preserved from erosion (see Fig. 1 for location) and sampled fluvial quartz-rich pebbles for ¹⁰Be profiles, which yielded an exposure age of 6400 ± 1300 yr. All details on dating are given in the GSA Data Repository item 2014152, available online (see data availability section).

We use the probabilistic approach developed by Zechar and Frankel²⁹ to compute slip rate distribution associated with terraces T2 (related to profiles P1) and T6 (profiles P2–P5 to the west and P6–P7 to the east). We assume normal distributions with a mean of 460 yr and 6400 yr and a
<figure><img src="image_15.png" /><figcaption>**Fig. 8** | Estimated trade-offs between model parameters for terrace T6. Same as Fig. 7, except we consider fault-perpendicular topographic profiles from P2 to P5.</figcaption></figure>

standard deviation of 70 yr and 1300 yr for T2 and T6, respectively, to model individual age uncertainties of these terraces. For cumulative displacement, we use two types of distribution¹: a boxcar to account for the range obtained with the rigid block model (Figs. 6 and 10b) and an arbitrary one provided by the pdf obtained by considering the thousand best-fitting models (Figs. 9 and 10c).

The computed slip rates are then compared to the convergence rate inferred from the present-day far-field interseismic velocity of 18 ± 2.1 mm.yr⁻¹ in central Bhutan³⁰. Regarding profile P1 across T2, the conventional model yields a slip rate lower than the convergence rate, which may be interpreted as a slip deficit and a likely higher seismic hazard (Fig. 10d). In parallel, our model suggests a slip rate higher than the convergence rate, which may imply a slip excess and a lower seismic hazard (Fig. 10e). However, the paleoseismic record available at the Sarpang site²⁰ indicates that T2 has probably only been affected by the 1714 CE earthquake since its deposition in the mid-16th century. This means that the time elapsed since deformation is an open interval, and any slip rate derived from it does not reflect the actual activity of the fault. Over the longer term, T6 has recorded at least five events in the last 2600 yr³¹ and likely 2–3 times as many since its deposition, which validates it as an adequate marker to characterize the Holocene deformation. When applied to profiles P2–P5 and P6–P7, the rigid-block model yielded slip rates of 18.2 ± 3.9 mm.yr⁻¹ and

20.2 ± 4.2 mm.yr⁻¹, respectively, and implied a significant variation over a distance of a few hundreds of meters (Fig. 10d). Conversely, our model is satisfied with a uniform slip rate of 19.6 ± 4.1 mm.yr⁻¹ (Fig. 10e). Although this variation may appear small with respect to uncertainties on slip rates (mostly due to high uncertainty on the age of T6), it remains that the rigid-block model would detect a lateral variation that would likely be interpreted as tectonic in origin (e.g. fault geometry, segmentation, slip vector, coupling), where our approach offers alternative explanations that should be explored (i.e., here a lateral variation in the thickness of the pre-deformation deposits).

## Conclusion

We propose an upgraded model to estimate offset from deformed fill terraces that underscores the main limitations of the conventional rigid block or kinematic approaches by accounting for abrupt changes in fault dip and curvature, mechanical properties of superficial material, and aggradation. We apply this model to the alluvial terraces observed in the Sarpang area, a region located in south-central Bhutan. The detailed mapping of alluvial terraces through field observations and analysis of a very high-resolution digital surface elevation model reveals the unevenness of terrace treads and a significant rugosity of the active river bed, with bars swelling reaching ca. 4 m elevation above the present stream level. In addition, geophysical
<figure><img src="image_16.png" /><figcaption>**Fig. 9 |** Model parameters and associated probability density function. Probability density functions (pdf) obtained from the 1000 best-fitted models out of the 2,000,000 stochastically explored. Red, green, and blue lines are related to fault-perpendicular topographic profiles P1, P2–P5 and P6–P7, respectively, from west to east. The fault geometry is defined from four dip angles α, β, γ, and δ, and lengths L₁, L₂, L₃, and L₄. The deepest segment is semi-infinite with L₄ = 100 km to avoid potential edge effects. The amount of slip *U*<sub>cumulative</sub> is constant along the fault. The</figcaption></figure>

slip mode of the fault, whether terraces are built with one or more earthquakes or by
creep is defined by the number of seismic events N. The pre-deformed topography
combines the current river slope of 0.9°, and the height H of the alluvial deposit. The
x-axis intervals for each parameter are based on the estimated range from geological
observations and geophysical investigations (see Supplementary Table 2 for more
details).

imaging suggests verticalization of the TFT beyond a depth of 5 m. Based on
this new information, modeling the highest terrace tread yields a slip rate of
19.6 ± 4.1 mm.yr⁻¹. The consistency of this value with the shortening rate
measured across central Bhutan confirms that most of the strain due to the
India-Eurasia collision is accommodated along the frontal thrust. In con-
trast, the fault slip rate obtained for the lowest terrace does not match the
convergence rate, suggesting that this terrace, deformed in the mid-16th
century, is not representative of Holocene deformation. These results
reassess the slip rate over two different periods and allow for a reevaluation
of seismic hazards in this highly populated region.

Methods

Terrace mapping and field survey

Our identification and characterization of fluvial terraces in the Sarpang
region is based on the combined analysis of remote sensing data and field
observations. Available digital elevation datasets for the Himalayas are
scarce, especially at high spatial resolution. Preliminary tests performed
on global datasets (e.g., SRTM, ASTER GDEM, AW3D30) publicly
available at 30 m resolution show noisy and ambiguous signals. In that
regard, we elected to produce an in-house very-high-resolution digital
surface model (DSM) derived from tri-stereo Pleiades satellite images by
photogrammetry using the MicMac software<sup>32</sup>. The resulting DSM
is extracted at 1 m resolution and comprises numerous off-ground points
from low vegetation, tree canopy, and buildings that need to be filtered
out to approximate a ground-point-only elevation model (Digital Ele-
vation Model). Hence, we applied a Cloth Simulation Filter (CSF)<sup>33</sup> as
implemented within CloudCompare software<sup>34</sup> with a multi-pass
approach to remove noise with decreasing frequency. CSF parameters
were especially fine-tuned to select and remove high vegetation observed
in the field (e.g., trees with heights ranging from 5 m to 20 m) while

preserving scarps. Furthermore, we applied an edge-preserving denoising
algorithm developed by Sun et al.<sup>35</sup>) to retain breaks-in-slope that may
help delimitate terraces. Where needed, we extracted topographic profiles
along 10-m-wide swaths by retaining the lowest elevation sampled
width-wise in order to filter out low vegetation (i.e., low trees and shrubs)
more efficiently. The final profiles (see Fig. 1a and in the “Result” section)
display few high-frequency low-amplitude pits and spikes (<2 m eleva-
tion) while retaining the low-frequency signal.

Fluvial terraces are defined here as abandoned surfaces corresponding
to the geomorphological markers formed by the ground surface, regardless
of their cut or fill nature. The initial identification is based on the analysis of
the digital elevation model (DEM) with the delimitation of continuous
surfaces exhibiting low slope (<2°), near-null curvature and low roughness
(see ref. 36) for a formalized approach) bounded by steep scarps from
hillslopes, erosion risers or faults. These surfaces were then visually con-
fronted with their local geomorphological context to ascertain they were
genetically connected with alluvial processes (e.g., proximity to drainage,
increasing median elevation with distance to drainage). Overall, we iden-
tified a total of eight alluvial terraces at the Sarpang site (see T0 to T7 in
Fig. 1a). They crop out at heights ranging from a few meters to ca. 60 m
above the mean elevation of the present-day stream water level (Fig. 1b). The
continuity of alluvial deposits suggests the terraces are nested. The height
distribution indicates they are unpaired, which may suggest rapid
downcutting.

Electrical resistivity tomography

Electrical Resistivity Tomography (ERT) is a geophysical method for im-
ging near-surface structures geometry (e.g. refs. 37–39). For this study,
electrical resistivity profiles are recorded along the eastern bank of the
Sarpang River (Fig. 1) with Wenner-Schlumberger (WS) and
<figure><img src="image_17.png" /><figcaption>**Fig. 10** | Sarpang terraces age, cumulative displacement along the Topographic Frontal Thrust, and associated slip rate distribution. a T2 and T6 ages distribution based on the Gaussian uncertainty model using data from refs. 12,20. b Displacement distribution derived from the conventional rigid block model, and a boxcar uncertainty model. Red, green, and blue lines correspond to the distribution obtained from profiles P1, P2–P5, and P6–P7, respectively. c Same as b using the</figcaption></figure>

proposed model based on an elastic dislocation full-space approach. d Slip rate distribution estimated from age and displacement modeled by a rigid block approach. The gray histogram shows the present-day far-field convergence rate distribution inferred from geodetic measurements in central Bhutan³⁰. e Same as d using our upgraded approach.

Dipole–Dipole (DD) configurations, and three different electrode spacings (1 m, 2.5 m, and 5 m) to take advantage of various investigation scales and resolutions²¹, proposed a resistivity model of the subsurface for each dataset with its own configuration and electrode spacing. However, this approach suffers from poor knowledge of the reliability of the obtained image, which is not only related to the electrode spacing and profile length but also depends on the distribution of the electrical resistivity at depth²³. In addition, this approach does not take full advantage of the dataset, particularly by combining high resolution and a large area of investigation. To overcome these issues, we develop a nested approach, where we invert the different surveys by integrating all the datasets in a single final acquisition file and following a three-stage inversion scheme. The inversion procedure starts with the most extended data profiles with 5 m electrode spacing to reconstruct the background resistivity model and large-scale structures. This model is then

refined by incorporating the intermediate spacing data, and finally, the complete acquisition file is inverted. The number of electrodes increases from 48 to 432 with denser measurements and deeper investigation depth around the fault zone. The PyGIMLi²⁴ and PyMERRY²³ codes are used for inversion and reliability assessment.

Regularization is commonly applied to avoid ill-conditioned inverse problems, but it causes smoothing of the ERT image, limiting the ability to define fault geometry accurately. To overcome this problem, we use a stochastic approach to improve accuracy following⁴⁰. We consider prior information from the ERT image and surface observation to set the range of parameters to test. We assume a simplified model that comprises two lithological units on either side of the fault. The fault dips northwards, and the dip angle can vary with depth. At the surface, we fix the geometry of the shallowest segment (fault location, dip angle, and length),
considering field observations. The model encompasses nine parameters, including the resistivity of each layer, the thickness of the two subsurface units, and three possible dip angle changes in depth (see Supplementary Table 1). For each set of parameters, we perform forward modeling with PyGIMLi, and the suitability of the related model is then estimated with an L2-norm misfit function from the collection of 100,000 sampled models. We finally obtain the posterior probability of the four-bodies resistivity model parameters (see Supplementary Fig. 1). Since we do not apply smoothing, the resulting posterior probability distributions provide an accurate and well-constrained description of the fault geometry at depth, including dip angles and segment lengths.

## Displacement model and inversion approach

In Bhutan, the only available estimate of the Holocene slip rate along the frontal thrust is based on uplift measurements of fluvial terraces¹² and structural observations of dip bedding in the Sarpang area²²,¹², converted the uplift rate $v_z$ into slip rate $v$ assuming a rigid block model (Fig. 2a), for which

$$v = \frac{v_z}{\sin(\alpha)^j} \quad (1)$$

where $\alpha$ is the surface fault dip. This conventional approach is useful but suffers from several limitations and approximations. First, the preserved terraces are viewed as planar landforms uplifted by sliding along a single-plane fault with a constant dip angle. Hence, this approach is not suitable for studying deformed terraces. Kinematic models partly overcome this limitation by proposing analytical formulations to study folded terraces associated with listric or curviplanar thrust fault models with a regular increase in dip toward the surface¹⁶. Second, in rigid block and kinematics models, the mechanical behavior of the medium is not considered, and fault geometry solely controls the surface deformation¹¹, propose an alternative approach in which surface deformation is the result of repeated earthquakes along a fault embedded in an elastic half-space. Third, in these previous approaches, the amount of uplift is estimated by comparing the altitude of the top of the fill terrace with the current profile of the river bed, assuming elevation variations solely along the downstream direction⁷,⁴¹. For alluvial rivers, the spatial and temporal interactions between flow rate and sediment load lead to complex morphodynamic processes shaping current river beds¹². The presence of natural levees, point bars, banks, or pools can lead to elevation variation up to several meters along transversal profiles⁴³. Neglecting to consider the topography shaped by depositional processes results in overestimating the river terrace uplift and the associated tectonic slip.

In our study, we compare the results from two end-member models: the conventional rigid block model used by Berthet et al.¹² and an upgraded model that we propose to replicate the deformation of alluvial terraces, as shown in Fig. 2b. Our model considers fault dip variations at depth and elastic properties. Due to unconsolidated material and a free surface, the resulting fault geometry near the surface can be more complex than a listric fault. Our model extends the previous kinematics approaches¹⁶ by accounting for abrupt fault dip and curvature changes. In addition, our approach focuses on forming terraces close to the fault. The correct consideration of topography is, therefore, crucial. Elastic half-space models⁴¹ do not allow for the incremental slip controlling progressive evolution of the distance between the fault and the top of the terraces, nor do they consider spatial variations in topography. Our approach extends the applicability of elastic models by allowing for initial topography related to sediment aggradation and the progressive evolution of the terraces' tread elevation.

We relate surface displacement to the underlying fault geometry using CuTDE, a Python translation and optimization of the original Matlab code developed by Nikkhoo and Walter⁴⁴. This approach uses triangular dislocation elements embedded in a homogeneous elastic full-space, including the effect of topography. We calculate displacements using a 2.5D approach. Four north-dipping segments give the fault geometry, each with its dip angle

and length (Fig. 2). The deepest segment has a semi-infinite length of 100 km to replicate the crustal-scale fault geometry⁴⁵ and avoid potential edge effects along kilometer-long profiles. The slip mode of the fault, whether terraces are built with one or more earthquakes or by creep, is defined by the number of seismic events $N$ to simulate the incremental slip controlling the progressive terrace's deformation. For the sake of simplicity, we assume a similar slip along the four segments of the thrust fault. The initial topographic profile is set by shifting the current river profile with a uniform factor $H$ associated with the topography inherited from river sediment deposits. Finally, we add the displacement from each seismic or creep event to calculate the cumulative displacement $U_{cumulative}$ and the associated tread elevation, considering both vertical uplift and horizontal tectonic advection. We neglect diffusion since the variations in slope on the tread surface are minimal and we are not examining incision and its impact on the data profiles.

We perform the inversion by stochastically exploring these model parameters in ranges defined by field observations and geophysical imagery (Supplementary Table 2). We test two million scenarios, for which we calculate their likelihood $L$ relative to the observations, defined as Eq. (2):

$$L = \exp\left(-\frac{1}{n}\sum_{i=1}^{n}(z_{i}^{\text{obs}} - z_{i}^{\text{calc}})^{2}\right) \quad (2)$$

with $n$ the amount of data along the studied profile, $z_{i}^{\text{obs}}$ the observed elevation and $z_{i}^{\text{calc}}$ the calculated one. The likelihood ranges from 0 (poor agreement with the observations) to 1 (perfect fit).

## Data availability

Elevation data is derived from Pleiades images provided under license by the French Space Agency (CNES). Topographic profiles used in this study are available in the following repository: https://doi.org/10.5281/zenodo.13851067. Resistivity data are available in the same publicly accessible repository https://doi.org/10.5281/zenodo.13851067. Geodetic data are already published³⁰ and given in table S1- supporting information (https://agupubs.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1002%2F2016GL071163&file=grl55338-sup-0001-supplementary.docx). Dating data are available in the paper of¹². Details are given in supplemental material (https://ndownloader.f庄static.com/files/23297183). Contact the corresponding author of this paper for resistivity data.

## Code availability

PyGIMLi²⁴ is an open-source Python library for multi-method modeling and inversion in geophysics. See the installation section to get this code (https://www.pygimli.org/installation.html). PyMERRY²³ is a Python code for post-inversion reliability assessment of electrical resistivity tomography images (https://doi.org/10.1190/geo2023-0105.1). CuTDE is a Python translation and optimization done by Ben Thomson from the Matlab code developed by Nikkhoo and Walter⁴⁴. This code is available through conda-forge (conda install -c conda-forge cutde) and can be found in the public github repository (https://github.com/tbenthompson/cutde). The MATLAB code used to produce the age, slip and slip rate distributions is based on the approach developed by Zechar and Frankel²⁹. It can be found in the auxiliary material of their paper (https://agupubs.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1029%2F2009JB006325).

Received: 10 December 2023; Accepted: 2 October 2024; Published online: 14 October 2024

## References

1. Ma, Z. et al. Tectonic and climate controls on river terrace formation on the northeastern Tibetan Plateau: Evidence from a terrace record of the Huangshui River. *Quat. Int.* **656**, 16–25 (2023).

2. Bookhagen, B., Fleitmann, D., Nishiizumi, K., Strecker, M. R. & Thiede, R. C. Holocene monsoonal dynamics and fluvial terrace formation in the northwest Himalaya, India. *Geology* **34**, 601–604 (2006).
3. Hassenruck-Gudipati, H. J., Ellis, T. S., Goudge, T. A. & Mohrig, D. A multi-proxy assessment of terrace formation in the lower Trinity River valley, Texas. *Earth Surf. Dyn. Discuss.* **2021**, 1–22 (2021).

4. Baynes, E. R., Lague, D. & Kermarrec, J. J. Supercritical river terraces generated by hydraulic and geomorphic interactions. *Geology* **46**, 499–502 (2018).

5. Hancock, G. S. & Anderson, R. S. Numerical modeling of fluvial strath-terrace formation in response to oscillating climate. *Geol. Soc. Am. Bull.* **114**, 1131–1142 (2002).

6. Molnar, P. et al. Quaternary climate change and the formation of river terraces across growing anticlines on the north flank of the Tien Shan, China. *J. Geol.* **102**, 583–602 (1994).

7. Lavé, J. & Avouac, J. P. Active folding of fluvial terraces across the Siwaliks Hills, Himalayas of central Nepal. *J. Geophys. Res. Solid Earth* **105**, 5735–5770 (2000).

8. Thompson, S. C. et al. Late Quaternary slip rates across the central Tien Shan, Kyrgyzstan, central Asia. *J. Geophys. Res. Solid Earth*, **107**, https://doi.org/10.1029/2001JB000596 (2002).

9. Simoes, M., Avouac, J. P. & Chen, Y. G. Slip rates on the Chelungpu and Chushiang thrust faults inferred from a deformed strath terrace along the Dungpuna river, west central Taiwan. *J. Geophys. Research. Solid Earth*, **112**, https://doi.org/10.1029/2005JB004200 (2007).

10. Burgess, W. P., Yin, A., Dubey, C. S., Shen, Z. K. & Kelty, T. K. Holocene shortening across the Main Frontal Thrust zone in the eastern Himalaya. *Earth Planet. Sci. Lett.* **357**, 152–167 (2012).

11. Hetzel, R., Niedermann, S., Tao, M., Kubik, P. W. & Strecker, M. R. Climatic versus tectonic control on river incision at the margin of NE Tibet: $^{10}$Be exposure dating of river terraces at the mountain front of the Qilian Shan. *J. Geophys. Res. Earth Surface*, **111** https://doi.org/10.1029/2005JF000352 (2006).

12. Berthet, T. et al. Active tectonics of the eastern Himalaya: New constraints from the first tectonic geomorphology study in southern Bhutan. *Geology* **42**, 427–430 (2014).

13. Bender, A. M. et al. Differential uplift and incision of the Yakima River terraces, central Washington State. *J. Geophys. Res. Solid Earth* **121**, 365–384 (2016).

14. Li, H. et al. Block tectonics across western Tibet and multi-millennial recurrence of great earthquakes on the Karakax fault. *J. Geophys. Res. Solid Earth* **126**, e2021JB022033 (2021).

15. Suppe, J. & Medwedeff, D. A. Geometry and kinematics of fault propagation folding. *Eclogae Geol. Helv.* **83**, 409–454 (1990).

16. Amos, C. B., Burbank, D. W., Nobes, D. C. & Read, S. A. L. Geomorphic constraints on listric thrust faulting: Implications for active deformation in the Mackenzie Basin, South Island, New Zealand. *J. Geophys. Res.* **112**, B03S11 (2007).

17. Charreau, J., Avouac, J. P., Chen, Y., Dominguez, S. & Gilder, S. Miocene to present kinematics of fault-bend folding across the Huerguosi anticline, northern Tianshan (China), derived from structural, seismic, and magnetostratigraphic data. *Geology* **36**, 871–874 (2008).

18. Coddington, J. A. & Burgette, R. J. Fold kinematics and fault slip rates from progressively deformed terraces: Implications for structural evolution of basins in the Kyrgyz Tien Shan. *Tectonics* **39**, e2019TC005776 (2020).

19. Trexler, C. C., Cowgill, E., Spencer, J. Q. & Godoladze, T. Rate of active shortening across the southern thrust front of the Greater Caucasus in western Georgia from kinematic modeling of folded river terraces above a listric thrust. *Earth Planet. Sci. Lett.* **544**, 1163–62 (2020).

20. Le Roux-Mallouf, R. et al. First paleoseismic evidence for great surface-rupturing earthquakes in the Bhutan Himalayas. *J. Geophys. Res. Solid Earth* **121**, 7271–7283 (2016).

21. Drukpa, D., Gautier, S., Cattin, R., Namgay, K. & Le Moigne, N. Impact of near-surface fault geometry on secular slip rate assessment derived from uplifted river terraces: Implications for convergence accommodation across the frontal thrust in southern Central Bhutan. *Geophys. J. Int.* **212**, 1315–1330 (2018).

22. Long, S., McQuarrie, N., Tobgay, T., Grujic, D. & Hollister, L. Geological map of Bhutan. *J. Maps* **7**, 184–192 (2011).

23. Gautier, M., Gautier, S. & Cattin, R. PyMERRY: A Python solution for an improved interpretation of electrical resistivity tomography images. *Geophysics* **89**, 1 (2024).

24. Rücker, C., Günther, T. & Wagner, F. M. PyGIMLi: An open-source library for modelling and inversion in geophysics. *Comput. Geosci.* **109**, 106–123 (2017).

25. Regmi, A. D. et al. Effect of rock weathering, clay mineralogy, and geological structures in the formation of large landslide, a case study from Dumre Besei landslide, Lesser Himalaya Nepal. *Landslides* **10**, 1–13 (2013).

26. Samouëlian, A., Cousin, I., Tabbagh, A., Bruand, A. & Richard, G. Electrical resistivity survey in soil science: a review. *Soil Tillage Res.* **83**, 173–193 (2005).

27. Giao, P. H., Chung, S. G., Kim, D. Y. & Tanaka, H. Electric imaging and laboratory resistivity testing for geotechnical investigation of Pusan clay deposits. *J. Appl. Geophys.* **52**, 157–175 (2003).

28. Holtmann, R., Cattin, R., Simoes, M. & Steer, P. Revealing the hidden signature of fault slip history in the morphology of degrading scarps, *Nature. Sci. Rep.* **13**, 3856 (2023).

29. Zechar, J. D., & Frankel, K. L. Incorporating and reporting uncertainties in fault slip rates. *J. Geophys. Res. Solid Earth*, **114** https://doi.org/10.1029/2009JB006325 (2009).

30. Marechal, A. et al. Evidence of interseismic coupling variations along the Bhutan Himalayan arc from new GPS data. *Geophys. Res. Lett.* **43**, 12–399 (2016).

31. Le Roux-Mallouf, R. et al. A 2600-yr-long paleoseismic record for the Himalayan Main Frontal Thrust (Western Bhutan). *Solid Earth* **11**, 2359–2375 (2020).

32. Rupnik, E., Daakir, M. & Pierrot Deseilligny, M. MicMac – a free, open-source solution for photogrammetry. *Open Geospatial Data Softw. Standards*, 2, https://doi.org/10.1186/s40965-017-0027-2 (2017).

33. Zhang, W. et al. An Easy-to-Use Airborne LiDAR Data Filtering Method Based on Cloth Simulation. *Remote Sens.* **8**, 501 (2016).

34. Girardeau-Montaut, D. CloudCompare (version 2.14) [GPL Software]. Retrieved from http://www.cloudcompare.org/ (2022).

35. Sun, X., Rosin, P. L., Martin, R. R. & Langbein, F. C. Fast and Effective Feature-Preserving Mesh Denoising. *IEEE Trans. Vis. Comput. Graph.* **13**, 925–938 (2007).

36. Iacobucci, G., Piacentini, D. & Troiani, F. Enhancing the Identification and Mapping of Fluvial Terraces Combining Geomorphological Field Survey with Land-Surface Quantitative Analysis. *Geosciences* **12**, 425 (2022).

37. Demanet, D. et al. The use of geophysical prospecting for imaging active faults in the Roer Graben, Belgium. *Geophysics* **66**, 78–89 (2001).

38. Nguyen, F., Garambois, S., Jongmans, D., Pirard, E. & Loke, M. H. Image processing of 2D resistivity data for imaging faults. *J. Appl. Geophys.* **57**, 260–277 (2005).

39. Villani, F. et al. Shallow subsurface imaging of the Piano Di Pezza active normal fault (central Italy) by high-resolution refraction and electrical resistivity tomography coupled with time-domain electromagnetic data. *Geophys. J. Int.* **203**, 1482–1494 (2015).

40. Ramirez, A. L. et al. Stochastic inversion of electrical resistivity changes using a Markov Chain Monte Carlo approach. *J. Geophys. Res.* **110**, B02101 (2005).

41. Wang, Y. et al. Deducing crustal-scale reverse-fault geometry and slip distribution from folded river terraces, Qilian Shan, China. *Tectonics* **39**, e2019TC005901 (2020).

42. Redolfi, M., Carlin, M. & Tubino, M. The Impact of climate change on river alternate bars. *Geophys. Res. Lett.* **50**, e2022GL02072 (2023).
43. Bristow, C. S. & Best, J. L. Braided rivers: perspectives and problems. *Geol. Soc. Lond. Spec. Publ.* **75**, 1–11 (1993).

44. Nikkhoo, M. & Walter, T. R. Triangular dislocation: an analytical, artefact-free solution. *Geophys. J. Int.* **201**, 1119–1141 (2015).

45. Le Roux-Mallouf, R. et al. Evidence for a wide and gently dipping Main Himalayan Thrust in western Bhutan, *Geophys. Res. Lett.*, https://doi.org/10.1002/2015GL063767 (2015).

## Acknowledgements

This work was supported by CNES, focused on the Pleiades satellite. Authors acknowledge financial support from the Agence Nationale de la Recherche (France): project TopoExtreme (PI: R.C. - grant ANR-18-CE01-0017) and the French Ministry for Higher Education and Research (M.G.'s PhD fellowship). Authors also thank Dr. Ben Thomson for access to CuTDE code and the Department of Geology and Mines of Bhutan for continuous logistic support during field campaigns. We express our gratitude to Y. Shang, J. Aslin, and the four reviewers for their valuable comments, which significantly enhanced the initial version of this paper.

## Author contributions

R.C., M.F., and S.G. initiated this study. M.G., S.G. and R.C. developed the ERT approach. M.F. produced the photogrammetric DSM and processed elevation data. R.C. developed the fault slip modeling approach. M.F., R.C., S.G., R.L., and D.D. performed field measurements. All authors contributed to the writing of the paper and provided critical feedback.

## 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/s43247-024-01759-z.

**Correspondence** and requests for materials should be addressed to Rodolphe Cattin.

**Peer review information** *Communications Earth & Environment* thanks Zhongtai He and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Yuan Shang and Joe Aslin. A peer review file is available.

**Reprints and permissions information** is available at http://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) 2024