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

Nearly ubiquitous along fluvial valleys, river terraces are major geomorphological archives of climate and tectonic processes (e.g. ref. ¹). 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. ²), and digital elevation model analysis to map these landforms³. 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. ⁴), and numerical modeling to study their emplacement through climatic modulation (e.g. ref. ⁵).

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,7,8,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^{8,10,11,12,13,14}. 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,16,17,18,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. ¹⁶ 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, 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.

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 terraces of the Sarpang river.

Full size image

Here, we propose an upgraded model for interpreting deformed alluvial 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 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 reevaluate the accommodation of deformation in this region, particularly during the Holocene.

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.

Full size image

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.

Full size image

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 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^12,20,21. 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.

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 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.

Full size image

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 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 (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^12,20 (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^23,24 and² a stochastic approach to account for resistivity contrast on both sides of the TFT²¹.

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 photogrammetric 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., \(\alpha =20^{\circ} -30^{\circ}\) for the shallow section and \(\beta =50^{\circ} -80^{\circ}\) for the deeper section.

Full size image

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^26,27. 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 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^12,20. 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).

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.

Full size image

Estimated trade-offs between the geometry of the steepest fault segment (dip angles \(\gamma\) and \(\delta\)), 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 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.

Full size image

Same as Fig. 7, except we consider fault-perpendicular topographic profiles from P2 to P5.

Full size image

The cumulative slip can also be estimated from Eq. (1) presented in “Methods” section for the flat treads 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.

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}_{{{{\rm{cumulative}}}}}\), 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}_{{{{\rm{cumulative}}}}}\). 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²⁸. Similarly, the pdf associated with the length \({L}_{1}\) and \({L}_{2}\) and dip \(\alpha\) and \(\beta\) 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.

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 \(\alpha\), \(\beta\), \(\gamma\), and \(\delta\), and lengths \({L}_{1}\), \({L}_{2}\), \({L}_{3}\), and \({L}_{4}\). The deepest segment is semi-infinite with \({L}_{4}\) = 100 km to avoid potential edge effects. The amount of slip \({U}_{{{{\rm{cummulative}}}}}\) is constant along the fault. 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\). 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).

Full size image

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).

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 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.

Full size image

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 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-16^th 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^-1, 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^-1 (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).

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 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^-1. 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 contrast, the fault slip rate obtained for the lowest terrace does not match the convergence rate, suggesting that this terrace, deformed in the mid-16^th 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.

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 tools³². 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 Elevation Model). Hence, we applied a Cloth Simulation Filter (CSF)³³ as implemented within CloudCompare software³⁴ 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.³⁵) 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 elevation) 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. ³⁶) for a formalized approach) bounded by steep scarps from hillslopes, erosion risers or faults. These surfaces were then visually confronted 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 identified 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 imaging near-surface structures geometry (e.g. refs. ^37,38,39,). For this study, electrical resistivity profiles are recorded along the eastern bank of the Sarpang River (Fig. 1) with Wenner-Schlumberger (WS) and 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^22,12. converted the uplift rate \({v}_{z}\) into slip rate \(v\) assuming a rigid block model (Fig. 2a), for which

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^7,41. 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}_{{{{\rm{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):

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