# 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

Dowchu Drukpa, Stéphanie Gautier, Rodolphe Cattin, Kinley Namgay, Nicolas Le Moigne

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

## ► To cite this version:

Dowchu Drukpa, Stéphanie Gautier, Rodolphe Cattin, Kinley Namgay, Nicolas Le Moigne. 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. Geophysical Journal International, 2018, 212 (2), pp.1315-1330. 10.1093/gji/ggx478. hal-02131196

HAL Id: hal-02131196

https://hal.umontpellier.fr/hal-02131196

Submitted on 11 Jun 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

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

# 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

Dowchu Drukpa,<fnref n="1" />,<fnref n="2" /> Stephanie Gautier,<fnref n="1" /> Rodolphe Cattin,<fnref n="1" /> Kinley Namgay<fnref n="2" /> and Nicolas Le Moigne<fnref n="1" />

Accepted 2017 November 6. Received 2017 September 18; in original form 2017 April 26

## SUMMARY

Vertical velocities obtained from uplifted river terrace dating near mountain fronts are commonly converted into overthrusting slip rates assuming simple geometry of the fault at depth. However, the lack of information on the dip angle of these shallow structures can lead to misinterpretation in the accommodation of convergence, and thus to erroneous conclusions on the transfer of shortening to the emergent thrust faults. Here, to assess the impact of fault geometry, we focus on the eastern Himalayan region in the south Central Bhutan, where the topographic frontal thrust (TFT) has been already documented by GPS, paleoseismic, geomorphic and geological studies. This study is based on high-resolution near-surface geophysical investigations, including electrical resistivity, seismic and gravity measurements. Using a similar stochastic inversion approach for all data sets, new quantitative constraints on both fault geometry and petrophysical parameters are obtained to image shallow depths, in the upper *ca.* 80 m. The combined results from both surface observations and geophysical measurement provide a TFT geometry that is dipping northwards with a shallow angle at the top (0–5 m), steeply dipping in the middle (5–40 m) and flattening at deeper depths (>40 m). Together, our new constraints on the fault geometry allow us to estimate a minimum overthrusting slip rate of 10 ± 2 mm yr⁻<fnref n="1" />, which is only a part of the *ca.* 17 mm yr⁻<fnref n="1" /> GPS convergence. This suggests that, in the study area, significant deformation partitioning on several faults including TFT and the Main Boundary Thrust cannot be ruled out. More importantly, assuming constant slip rate, the obtained dip angle variations lead to uplift rate changes with distance to the TFT. This underlines that taking into account uplift rate from terrace dating only at the front location and assuming a constant dip angle fault geometry based on surface observations may significantly bias the slip rate estimates.

**Key words:** Geomorphology; Gravity anomalies and Earth structure; Electrical resistivity tomography (ERT); Seismic tomography.

## 1 INTRODUCTION

The Himalaya that stretches *ca.* 2500 km from the Hazara-Kashmir syntaxis in the west to the Namcha Barwa syntaxis in the east constitutes one of the most seismically active regions of the world. In this area, many previous studies (e.g. Molnar & Tapponnier 1975; Seeber & Armbruster 1981; Pandey *et al.* 1995; Bilham *et al.* 1997; Cattin & Avouac 2000; Lavé & Avouac 2000; Decelles *et al.* 2002; Kumar *et al.* 2010; Duputel *et al.* 2016) have shown that the occurrence of major earthquakes can be related to the India-Eurasia shortening, which has been accommodated along the Main Himalayan Thrust (MHT), a mid-crustal decollement where the India plate is underthrusted beneath the Himalayas and Tibet.

In central Nepal and Arunachal Pradesh, estimates of Holocene horizontal shortening rates have been already obtained from studies of uplifted river terraces (Lavé & Avouac 2000; Burgess *et al.* 2012) nearby the Main Frontal Thrust (MFT), which is the most recent surface expression of the MHT (Schelling & Arita 1991; Pandey *et al.* 1995). In these areas, it has been reported that the MFT absorbs most of the shortening rate across the Himalaya, whereas the MFT is locked over interseismic periods (Ader *et al.* 2012). This leads to the current understanding of seismic cycle in central Himalaya where most of the interseismic deformation deficit is released during *M* > 8 earthquakes that rupture the MFT up to the surface (Bilham *et al.* 1997; Cattin & Avouac 2000).

Located in the eastern part of the Himalayas, the Kingdom of Bhutan remains one of the least studied segments of the

<footnote n="1">¹Geosciences Montpellier, UMR5243, Université de Montpellier, Place E. Bataillon, F-34095 Montpellier, France. E-mail: cattin@gm.univ-montp2.fr</footnote>

<footnote n="2">²Seismology and Geophysics Division, Department of Geology and Mines, Ministry of Economic Affairs, Post Box 173, Thimphu, Bhutan</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Location of the Kingdom of Bhutan and the study area in south central part of Bhutan. (b) High-resolution Pleiades image of the study area showing the Topographic Frontal Thrust (TFT) fault trace, the location of the paleoseismic trench studied by Le Roux-Mallouf *et al.* (2016) and our geophysical profile (yellow line). (c) Seismic survey layout plan.</figcaption></figure>

Himalayan range in terms of active tectonics and its associated seismic risk assessment (Drukpa et al. 2006). Over the last decade several geological and geophysical studies have been carried out in Bhutan to fill this information gap (e.g. Hammer et al. 2013; Berthet et al. 2014; Vernant et al. 2014; Le Roux-Mallouf et al. 2015, 2016; Hetényi et al. 2016a,b; Marechal et al. 2016; Diehl et al. 2017; Singer et al. 2017). In particular, a multidisciplinary approach including geomorphology, paleoseismology and geodesy has been conducted in Sarpang, a small town in south Central Bhutan (Fig. 1). In this area, these studies suggest a seismic behaviour similar to the one observed in other parts of Himalayas, with (1) present-day deformations indicating a fully locked thrust-fault at the surface over the interseismic period (Marechal et al. 2016) and (2) evidences of occurrence of major earthquakes over the last millennium (Berthet et al. 2014; Le Roux-Mallouf et al. 2016).

Tectonic scarps and well-preserved abandoned terraces on both banks of the Sarpang river also attest to the accumulation of vertical deformation through time (Fig. 1). Based on both radiocarbon and cosmogenic ¹⁰Be dating, Berthet et al. (2014) estimate a Holocene vertical uplift rate of 8.8 ± 2.1 mm yr⁻¹ for these terraces. By assuming a constant dip angle of 25° ± 5° towards the north, Berthet et al. (2014) calculated a Holocene slip rate of ~20 mm yr⁻¹ along the main thrust fault in Central Bhutan, suggesting that the entire convergence rate across Bhutan Himalaya is accommodated at the front without any deformation partitioning. However the most recent

geological map of Bhutan proposed by Long *et al.* (2011) presents dip angle values in ranges of 10°–60° for the Sarpang area, making slip rate difficult to assess from fault geometry deduced from surface observations only.

Additional constraints on the geometry are thus required to estimate the amount of shortening rate accommodated along the frontal fault in Bhutan, which is a key input parameter for seismic hazard assessment. Here, we image the shallow structures from high-resolution near-surface geophysical methods including electrical resistivity, seismic refraction and gravity measurements. Both electrical and seismic methods are proved to be essential tools in active fault imaging by providing new constraints on fault geometry, fault offset and abrupt variations of sediment thickness (e.g. Morey & Schuster 1999; Storz et al. 2000; Stephenson & McBride 2003; Karastathis et al. 2007; Nguyen et al. 2007; Dorn et al. 2010; Villani et al. 2015). In comparison with the paleoseismic trench studied by Le Roux-Mallouf et al. (2016), these two methods allow for a wider investigation area and larger exploration depth. Because gravimetry takes into account large-scale anomalies (Benson & Mustoe 1995; Wise et al. 2003), we also conduct a gravity survey in order to investigate deeper structures.

After a brief presentation of the study area, we describe data acquisition. Next we present inversion approaches implemented to assess the frontal fault geometry. Finally we discuss our results and implications in terms of shortening accommodated along the frontal fault.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Picture of excavated site in Sarpang area with inset in the bottom left corner showing the TFT trace at the surface. A sketch of the outcrop shows the estimated dip angle of the fault as observed in the trench.</figcaption></figure>

## 2 STUDY AREA

The Kingdom of Bhutan, located in the eastern part of the Himalaya, accounts for approximately 14 per cent of the total length of the Himalayan arc. Our study area is approximately 90 km southeast of the capital city, Thimphu. It is located in the Himalayan foothills along the Sarpang river in south Central Bhutan, close to the border with India (Fig. 1).

While the MFT is well-defined in Nepal and coincides with the present-day Himalayan topographic front, its location in Bhutan is less straightforward to define (Le Roux-Mallouf et al. 2016). In the Sarpang area where a 10-km-wide re-entrant feature directly juxtaposes lesser Himalaya and Quaternary alluvium (Long et al. 2011), the most frontal structure expressed in the geomorphology is an anticline that likely accommodates a limited, though, undetermined fraction of the shortening (Dasgupta et al. 2013). For consistency with previous studies, the main topographic scarp studied here will be referred to as the Topographic Frontal Thrust (TFT) of Bhutan (Fig. 1).

This site has been selected for several reasons. First, in many parts of Bhutan the TFT is located at the Indian border. The study site is in the Sarpang re-entrant, where the topography front is 10 km farther north compared to other foothill areas (Long et al. 2011). Furthermore, since the site is located along a river in the foothill, it is characterized by low elevation variations. This specific geographic feature ensures easy field accessibility and feasibility. Second, a host of information is now available regarding geology, geomorphology and paleoseismology for this area. The main lithological units in the fault zone documented on the trench log (Le Roux-Mallouf et al. 2016) comprises highly deformed phyllite in the deepest exposed unit on the hanging wall that corresponds to the Paleozoic Buxa Formation of the Lesser Himalaya Sequence. An alluvial layer is overlying this unit. Southward, in the footwall side, the gravel unit observed in the hanging wall is underlain by a

fine-grained sand layer deposited by the lateral drainage system. A shear-zone juxtaposes the hanging wall and footwall along a 1 m deformed zone comprising sheared sand and gravel unit. In the trench, this sheared and highly deformed zone with a dip angle of 17°–24° is interpreted to be the TFT (Fig. 2). Abandoned fluvial terraces and radiocarbon dating suggest a Holocene uplift rate of 8.8 ± 2.1 mm yr⁻¹ (Berthet et al. 2014) in this area. As previously mentioned, the secular slip rate assessment is however not straightforward because near the front in Sarpang and vicinity area, bedding and foliation of rocks exhibit high variations in their dip angles, with values ranging from 10° to 60° (Long et al. 2011).

Together, these features make the Sarpang area a relevant site for geophysical investigations to assess the geometry of the TFT at shallow depth.

## 3 DATA

We acquired geophysical data along the east side of Sarpang river where the east-west trending TFT trace was intersected by the previously mentioned paleoseismic trench down to 1 m depth (Figs 1 and 2). Taking into account the geological setting and the depth of the target, geophysical investigations included 2-D electrical tomographies, seismic refraction tomographies, and gravimetry measurements. The use of different geophysical methods allow to image different physical properties at different scales of investigation. All geophysical data were collected along the same N-S profile, with different spread lengths depending on the methods. The midpoints of the geophysical surveys were positioned at the fault location at the surface deduced from the paleoseismic study. We observed little and smooth topographic variation (<1.3 m) along the geophysical line and therefore envisaged no topographic correction for the seismic and resistivity data analyses.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** ERT pseudo-sections (Wenner–Schlumberger array) observed along the study profile for electrode spacing of 1 m (top), 2.5 m (middle) and 5 m (bottom).</figcaption></figure>

## 3.1 Electrical resistivity tomography (ERT)

We acquired three ERT profiles using 48 electrodes with 1, 2.5 and 5 m electrode spacing in order to combine high resolution at shallow depth (1 m spacing) and large investigation depth (5 m spacing) using an Iris Syscal multielectrode system. We adopted the commonly used Wenner–Schlumberger and dipole–dipole configuration to achieve a good compromise between vertical and horizontal resolution and effect of noise (Dahlin & Zhou 2004; Loke 2015). The resistivity values obtained from the measurement were checked for any inconsistency and cleaned. Both the Wenner–Schlumberger and dipole–dipole arrays provided consistent results. Hereinafter, we present only Wenner–Schlumberger electrical images because of a greater sensitivity to both lateral and vertical variations (Nguyen et al. 2007). The obtained ERT pseudo-sections show high electrical resistivity contrasts (~1:100) across the fault zone with a nearly vertical contact down to ~40 m depth (Fig. 3). Vertical resistivity changes are also highlighted at shallow depth south of the fault, underlying a high sensitivity of the Quaternary layers to electrical properties. The fault zone is characterized by an apparent resistivity of 100–1000 Ω m. The north side shows a uniform apparent resistivity layering with a thin upper layer resistivity of 200–1000 Ω m overlying a layer with a very low resistivity <100 Ω m. The south side shows relatively constant apparent resistivity values (1000–

4000 Ω m) with a very high resistivity zone located at 5–15 m depth at the southern end of our profile. A thin upper layer of low resistivity is also observed southward.

## 3.2 Seismic tomography

A seismic survey coincident with the ERT profiles was carried out in order to obtain a near-surface velocity model around the fault zone. We used a MoHo® Soilspy Rosina digital multichannel seismic acquisition system with 24 4.5-Hz vertical component geophones to record *P*-waveform traces to deduce *P*-wave velocity profile. A 1 m receiver spacing and five roll-alongs (shift of 20 geophones and overlap of 4 geophones each time) were used to finally acquire a 103 m long seismic profile. We shot with a 10 kg sledgehammer every 1 m along each seismic line. We recorded inline shots at every geophone positions and also offset shots starting 25 m before the first trace and ending 20 m after the last geophone. An average of 72 shots was recorded along each line resulting in a total of 360 shots along the 103 m long profile (Fig. 1). Our data have good quality with clear first breaks on waveform traces. We hand-picked a total of 5760 first-arrival traveltimes (Fig. 4). First-arrival times show smooth changes at ~30 m north of the fault trace. The northern part of the profile presents small values of first arrival times (~10 ms)
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Hodochrone of the hand-picked first-arrival traveltimes associated with the roll-along seismic acquisition. Colour dots correspond to a first-arrival time for a given geophone positions along the profile. White dots correspond to the case for which shot and station are at the same location along the seismic profile. The colour scale represents the value of first-arrival times.</figcaption></figure>

whereas an abrupt transition between small (~10 ms) and higher (~25 ms) arrival times is observed in the south, indicating the presence of a shallow interface at depth on this side.

## 3.3 Microgravity

In addition to electrical resistivity and seismic investigations, we performed microgravity measurements using Scintrex CG-5 gravimeter (~10 µGal resolution) along the same north-south profile. From the centre point of the survey line positioned at the paleoseismic trench, we collected gravity readings every 5 m on the either side of the profile covering a distance of 30 m and 105 m to the south and to the north, respectively. Spentially denser gravity points every 1 m were also collected in the vicinity of the fault area. Using the GravProcess software (Cattin et al. 2015) network adjustment was performed and topographic effect was corrected from accurate elevation data gathered along the same profile assuming a constant density of 2670 kg m⁻³. A regional trend of -1.58 µGal m⁻¹ obtained by Hammer et al. (2013) is also taken into account. Finally our data set consists of 139 corrected gravity measurements, which highlight variations along the profile. No change at the fault trace is observed but a transition occurs at around 27 m north of the TFT. The southern part of the profile is characterized by a moderate northward increase of ca. 4 µGal m⁻¹. The northern part shows an increase twice as large with a northward increase of ca. 450 µGal in 65 m (Fig. 5).

# 4 INVERSION APPROACH

Following Ramirez *et al.* (2005), we use a stochastic method for imaging shallow structures from the obtained geophysical data sets.

The approach combines prior information and forward modelling to produce models consistent with the available data. Following Mosegaard & Tarantola (1995), a Markov Chain Monte Carlo technique is used to pseudo-randomly generate a large collection of models according to the posterior probability distribution.

Assuming a simplified geometry, each model is associated with only five bodies including a south (STL) and north (NTL) top layers, a south (SL) and north (NL) shallow layers and a fault layer (Fig. 6). On the basis of this formulation, a model combines 10 parameters, including the velocity, the resistivity and the density of each body as well as the thickness of layers, the fault location and the fault dip angle.

A priori parameter ranges for resistivity and velocity values are obtained from preliminary deterministic inversion using Res2dInv (Loke & Barker 1996) and Rayfract (Schuster & Quintus-Bosz 1993; Sheehan et al. 2005; Pasquet et al. 2015), respectively (Supporting Information Figs S6 and S7). We assume prior density contrasts between NL and the other bodies in -500 to 500 kg m⁻³ range Concerning the geometry prior information comes from structural and geomorphic observations (e.g. Long et al. 2011; Le Roux-Mallouf et al. 2016), which yield top layer thickness less than 5 m and a fault dip angle between 10° and 80° (see Supporting Information Table S1).

The pseudo-random walk through this multidimensional parameters space is controlled by the following rules for the transition between model $m_i$ to model $m_j$:

(1) If $L(m_j) \ge L(m_i)$ then accept the proposed transition from $i$ to $j$.

(2) If $L(m_j) < L(m_i)$ then accept the proposed transition from $i$ to $j$ with the probability $\frac{L(m_j)}{L(m_i)}$,
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Measured gravity (top), elevation (middle) and gravity variations corrected for both topographic effect and regional trend (bottom) along the study profile. Data uncertainty is associated with both accuracy of the CG5 gravimeter and error in elevation measurement.</figcaption></figure>

where $L(m_i)$ and $L(m_j)$ are the likelihood of the old and the new model, respectively. Here we assume that the likelihood function can be written as

$$L(m_i) = \exp \left( -\frac{1}{n_{\text{obs}}} \sum_{n=1}^{n_{\text{obs}}} \frac{|\text{calc}_n(m_i) - \text{obs}_n|}{\sigma_n} \right)$$

where $n_{\text{obs}}$ is the number of data points, ‘obs’ is the data vector and $\sigma^2$ is the total variance, that is, the uncertainties associated

with each data point. calc ($m_i$) is the forward modelling function associated with the model $m_i$. Depending of the considered data set, this function is obtained using the following forward formulations:

(1) 2-D geoelectrical modelling is performed with the software package R2 (Binley & Kemna 2005; Binley 2015). The current flow between electrodes is obtained using a quadrilateral mesh. The number of nodes in the vertical direction is fixed to 50 with an
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Geometry of the model used in the stochastic inversion. STL—South Top layer, NTL—North Top Layer, SL—South Layer and NL—North Layer. xfault is the prior location of the fault as observed in the field. Model thickness is associated with the thickness investigated by each geophysical method.</figcaption></figure>

exponentially increasing node spacing from 0.05 m at the surface to 45 m at 300 m depth. In the horizontal direction a constant node spacing of 0.5 m is used along the study profile.

(2) Synthetic traveltimes are computed using the real receiver-shot configuration and solving the Eikonal equation with a finite-difference algorithm (Podvin & Lecomte 1991). Rays are traced in the obtained time field with the *a posteriori* time-gradient method. More precise traveltimes are then estimated along ray paths (Gautier *et al.* 2006; Priollo *et al.* 2012). The velocity model is a 106 m × 5 m area and extends to 10 m depth in order to include all ray paths (Supporting Information Fig. S7). The model is discretized on a regular grid of 425 × 21 × 45 nodes with a spacing of 0.25 m and the velocity field is parametrized by trilinear interpolation between grid nodes.

(3) Gravity variations along the profile are calculated from the formulations developed by Won & Bevis (1987), which provide the gravitational acceleration due to *n*-sided polygons. Here the polygons are associated with the geometry of the five bodies described above. The model is extended to 10 km southwards and northwards to avoid edge effects at the two terminations.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Misfit between observed and calculated ERT pseudo-sections for electrode spacing of 1, 2.5 and 5 m using Wenner-Schlumberger configuration. This misfit is defined as the difference between the observed and calculated resistivities using the electrical model plotted at the bottom. Bottom colour scale is the same as Fig. 3 for comparison.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Distribution of TFT dip angle from ERT sections using both dipole-dipole and Wenner-Schlumberger arrays. Electrode spacing ranges from 1 m (top) to 5 m (bottom).</figcaption></figure>

The posterior probability of each model parameter is then obtained from the final collection of the 5.10⁵ sampled models (Supporting Information Figs S1–S5). Compared to commonly used approaches based on the search for the simplest model, the main advantages of our method include its ability (1) to assess the fault geometry because no smoothing is applied, (2) to provide a measurement of the uncertainties on the obtained dip angle and (3) to allow trade-off analysis between geometric and either electrical resistivity, velocity or density properties. Using parallelism, the computation time associated with electrical, seismic and gravity inversion on a 10 core work station is *ca.* 10 hr, 1 d and 30 min, respectively.

# 5 RESULTS

## 5.1 Electrical resistivity tomography

The set of most likely models derived from the stochastic approach explains the main features of the observed apparent resistivity pattern (Fig. 7) and points out a high fault dip angle of *ca.* 70° (Fig. 8). This is obtained for both dipole-dipole and Wenner-Schlumberger configuration as well as for electrode spacings between 1 and 5 m. Bivariate frequency histograms indicate no trade-off between dip

angle and the other geometric and electrical parameters (Supporting Information Figs S1–S3). Those histograms suggest a 2.5 m thick fault zone. However the resistivity of this unit remains poorly resolved. The inversion approach images thin low-resistive top layers, both on the southern (~2.5 m, ~550 Ω m) and northern sides (~3.5 m, ~350 Ω m) (see Supporting Information Table S1 for details). The small resistivity contrasts between those two top layers can prevent the estimation of the fault geometry at very shallow depth (<5 m). In contrary, due to the very high resistivity contrast between the two deeper bodies (SL ~ 3300 Ω m versus NL ~ 30 Ω m), we consider the obtained fault dip angle as a well-constrained parameter down to 40 m depth. This is confirmed by the narrow posterior distribution obtained for dip angle (Fig. 8). Finally, one can mention some discrepancies between the observed and calculated pseudo-sections on the south part (Fig. 7). Based on the simplified geometry of the model assuming horizontal layering, our inversion procedure is not able to explain the north-south resistivity variations in the footwall of the TFT.

All together these information on both geometry and resistivity contrast suggest a very clear resistivity contrast between the two sides of the fault as well as a constant dip angle of ~70° over a depth ranging between *ca.* 5 m and *ca.* 40 m.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Hodochrone of the residuals between observed data and theoretical first-arrival times computed in the velocity model plotted at the bottom. Compared to the resistivity model Fig. 7, note the lower limit of the y-axis, which corresponds to a shallower depth of investigation.</figcaption></figure>

## 5.2 Seismics

The set of final velocity models (Fig. 9) resulting from the stochastic inversion approach provides low traveltime residuals of ±3 ms for most of the source-receiver pairs. This suggests that our simple geometry captures the main features of the velocity field. Furthermore, traveltime residuals show a quite homogeneous pattern, except close to the fault trace between −5 and 15 m where residuals abruptly increase from −5 to 4 ms northward. This demonstrates that seismic data are influenced by the presence of the fault.

Ray coverage (Supporting Information Fig. S7) indicates a shallower resolution depth compared to ERT investigations. Resolution depth varies between the two sides of the fault, from *ca.* 8 m to *ca.* 5 m in south and north, respectively. At these depths, the velocity models resulting from both the stochastic inversion (Fig. 9) and tomography (Supporting Information Fig. S7) point out high velocity variations of *ca.* 50 per cent at the transition of the fault zone. Fig. 9 also emphasizes strong vertical velocity changes on both sides of the fault. High velocity contrasts between top and bottom layers induce a concentration of rays (Supporting Information Fig. S7) at

a depth between 2 and 4 m, which prevents deeper investigations, in particular in the north.

Taking into account this shallow investigation depth, the velocity field can be characterized by two deeper units of *Vp* ~ 1100 m s⁻¹ (SL) and *Vp* ~ 2100 m s⁻¹ (NL) below two superficial low-velocity layers (STL: ~5 m, ~900 m s⁻¹ and NTL: ~3 m, ~1600 m s⁻¹) (Fig. 9). Our inversion procedure also reveals that seismic data are sensitive to the dip angle parameter (Fig. 10). Our result suggests a northward dipping fault with a low-angle of *ca.* 20°–30° at depths down to *ca.* 5 m, which is consistent with field observations in the trench (Le Roux-Mallouf *et al.* 2016). This observation is also in agreement with both pseudo-sections and ERT profile (Fig. 3 and Supporting Information Fig. S6), which displays a change in dip angle with a more gentle slope of the TFT fault near the surface. Because we assumed a constant fault dip angle and we used the same model for the different ERT configurations, the inversion procedure was not able to image this dip angle change near the surface with resistivity data only.

Hence, field observations with the obtained seismic and electrical resistivity images together suggest a northward dipping fault with a
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Distribution of TFT dip angle obtained from seismic measurements.</figcaption></figure>

low angle of *ca.* 20°–30° at very shallow depth, increasing gradually up to *ca.* 70° at a depth of 5–10 m.

## 5.3 Gravity

The observed short wavelength decrease at 27 m (Fig. 11) from the fault is related to topography correction. As previously mentioned, gravity measurements are mostly affected by the deeper part of our model. Thus in the following we focus on the long wavelength of the gravity profile associated with the two south-north gravity gradients described in Section 3.3.

The result of our inversion suggests that the observed northward increase of gravity measurements is mostly related to both $\Delta\rho$, the density contrast between SL and NL, and $\alpha$, the fault dip angle (Figs 11 and 12). As indicated in Supporting Information Fig. S5, gravity measurements cannot be used to assess the other density and geometric parameters, which remain poorly constrained.

Our result reveals a clear trade-off between $\Delta\rho$ and $\alpha$: the higher the density contrast, the lower the fault dip angle. For $\Delta\rho = -350 \text{ kg m}^{-3}$ the fault dip angle is *ca.* 30°, whereas for $\Delta\rho = -200 \text{ kg m}^{-3}$ the fault dip angle is *ca.* 60° (Fig. 11). This leads to a wide distribution of fault dip angle (Fig. 12). The maximum of the distribution obtained at $\alpha \sim 30°–40°$ and for a model thickness ranging between 70 and 90 m suggest however a fault that flattens at depths below the investigation depth of 40 m obtained with ERT.

## 6 DISCUSSION

The geometry of the fault, especially at shallow depth, is a key parameter in estimating the slip rate on the fault from uplifted terrace dating. Here we discuss our account on the new constraints for the TFT geometry deduced from near-surface geophysical techniques

and its implication on convergence partitioning at the frontal thrust zone in south Central Bhutan.

### 6.1 Subsurface imaging

We take advantage of the various scales of investigation coming from ERT, seismic and gravity methods to obtain an accurate description of shallow structures and fault geometry at depth, which can be subdivided in three main zones: (1) a very shallow part up to 5 m depth well-constrained by both field observations and seismic data considering the ray coverage; (2) an intermediate depth part well-imaged by ERT sections between 5 and 40 m depth due to high-resistivity contrasts; (3) a deeper part documented by gravity measurements below 40 m depth. The fault geometry discussed here arises from the integration of these three scales.

The obtained fault thickness remains poorly constrained. Assuming an *a priori* thickness between 0 and 5 m, the gravity approach gives no information on this parameter and both the electrical and the seismic studies suggest a fault thickness of 2–5 m. In terms of lithological setting and water content, the geophysical data sets suggest a thin layer (~3–5 m) that appears to be present on both sides of the fault trace and which probably corresponds to recent alluvial deposits. Along the profile, resistivity and velocity variations at shallow depth may probably be due to a northward decrease of water saturation. Below these superficial layers, in the hanging wall of the TFT, the obtained very low-resistivity values of <30 'Ω m, the high $Vp$ of *ca.* 2100 m s$^{-1}$ and the relatively low densities can be associated with a phyllite unit, which can be observed in the field. Overall, the geophysical methods image a more complex fault geometry than proposed by earlier studies (Berthet *et al.* 2014; Le Roux-Mallouf *et al.* 2016). Our results show a TFT with a flat and listric-ramp geometry with a low dip angle of 20°–30° at shallow depth, steeply dipping at ~70° in the middle and gradually flattening in a shallower dip angle of 30°–40° in its deeper part (Fig. 13).
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top: comparison between observed (blue circles) and calculated (grey lines) gravity variations along the study profile obtained for the 100 best-fitting models. Bottom: density contrast models associated with the red and green lines plotted above and corresponding to fault dip angles of 30° and 60°, respectively.</figcaption></figure>

## 6.2 Overthrusting slip rate assessment

Berthet et al. (2014) estimate a Holocene vertical slip rate of $8.8 \pm 2.1 \text{ mm yr}^{-1}$ by dating two uplifted river terraces in the Sarpang area. Assuming a dip angle of $20°-30°$, they propose a slip rate of $20.8 \pm 8.8 \text{ mm yr}^{-1}$, which is consistent with the GPS convergence rate of 17 mm yr$^{-1}$ obtained across central Bhutan (Marechal et al. 2016). Finally they conclude that the Himalayan convergence is mostly accommodated by the TFT. However, this major conclusion can be revisited from our new constraints on the TFT geometry.

First, assuming a constant overthrusting slip rate along the TFT, a vertical velocity profile is calculated from this observed uplift rate (Okada 1985). As expected this calculated profile depends on TFT geometry (Fig. 14a): a higher fault dip angle implies a higher uplift rate. More surprisingly, it also depends on the distance between the TFT and the location of dating samples. For instance, a distance of 5 m with respect to the TFT yields two very different

vertical velocity profiles associated with the two end-member models for the fault geometry (Fig. 14a). On the contrary, if uplift rate is measured about 10 m north from the TFT, the uplift rate difference drastically reduces (Fig. 14b). In other words, due to the flat and listric-ramp geometry of the shallow TFT, the uplift rate measured on the top of river terraces are spatially variable and cannot be constant. This result questions the validity of commonly used approaches for which a mean uplift rate is obtained by combining several uplifted terraces located at various distances from the front. Furthermore, assuming that the far-field GPS shortening rate corresponds to an upper limit for the uplift rate, our calculation shows that part of the models are unrealistic (see green curve on Fig. 14a). This suggests that both the convergence rates derived from GPS and the uplift measurements can be used to reduce the *a priori* geometric parameter ranges tested in our stochastic approach.

Second, assuming no prior information on the relative location of uplift rate measurements, one can deduce the overthrusting slip
<figure><img src="image_15.png" /><figcaption>Figure 12. Distribution of TFT dip angle obtained from gravity measurements.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13. Simplified cross-section showing the final model obtained from electrical resistivity tomography, seismic refraction and gravity measurements. Together geophysical results suggest a TFT with a flat and listric-ramp geometry, with a surface dip angle of ~20° reaching ~70° at 20 m depth and flattening in its deeper part. Note that the dashed area is bounded by the two end-member models of fault geometry given by the green and the red lines. Hence, this area does not represent the fault thickness, which is estimated to be 2.5 m.</figcaption></figure>

rate from the TFT geometry (Fig 14c). The slip rate associated with a rigid block model with a constant dip angle α can be easily estimated from:

$$slip\ rate = \frac{uplift\ rate}{\sin(\alpha)}$$

As previously proposed by Berthet et al. (2014), this simple approach gives a minimum dip angle of 30° for which most of the convergence across central Bhutan is accommodated along the TFT (Fig. 14c). However, the steeper is the dip angle, the greater is the chance for slip partitioning with other faults. Assuming a constant uplift rate of 8.8 ± 2.1 mm yr⁻¹ along the profile, corresponding to
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Estimated overthrusting slip rate inferred from both fault geometry and observed uplift rate. (a) Uplift rate along the study profile. Red and green curves are associated with the two end-member models obtained for fault geometry described Fig. 13. Thick blue line denotes the far-field shortening rate estimated from GPS measurements (Marechal *et al.* 2016). It corresponds to the upper limit of uplift rate, which can be associated with a theoretical vertical fault. Grey circle is the observed uplift rate assuming a northward distance of 5 m from the TFT as reported by Berthet *et al.* (2014). (b) Same as (a) except a northward distance of 10 m is assumed. (c) Estimated slip along the TFT at depth assuming a constant uplift rate along the study profile for the two-end member models denoted by red and green curves. Hatched area around these curves is associated with uplift rate uncertainties. The thick grey dashed lines point out the area of uplift rate assuming a northward distance from the TFT of 5 and 10 m, respectively. Note that within this area, the uncertainties in the uplift rate spikes up close to the TFT and decreases away from the front towards the north. Thick blue line denotes the far-field shortening rate estimated from GPS measurements (Marechal *et al.* 2016), which is the upper limit of slip rate. The slip obtained from rigid blocks model with a constant dip angle ranging from 10° to 60° is given by black dashed lines. (d) Fault geometry used.</figcaption></figure>

the hypothesis of no information on the sampling location, the slip rate can also be estimated from a less straightforward modelling based on the obtained geometry. In that case using dislocations embedded in a homogeneous half space (Okada 1985), the obtained slip rate exhibits high variations along the profile from 20–

40 mm yr⁻¹ above the very shallow part of the fault to 10–20 mm yr⁻¹ in the northern part of the profile. Using the convergence rate as a maximum value for the slip, this result suggests a minimum distance of 8 m for the steepening of the TFT and an accommodation of at least 10 ± 2 mm yr⁻¹ of the 17 mm yr⁻¹ of convergence
on the TFT. The obtained uncertainties associated with this slip rate estimate arise mainly from the location of samples for terrace dating and the fault geometry inferred from geophysical inversion.

## 6.3 Deformation at the topographic front

Based on our new constraints on the TFT geometry and the resulting slip rate, we propose that at least 60 per cent of the convergence rate due to ongoing underthrusting of India beneath the Himalaya is accommodated by the TFT. It results that additional faults must be active in this area, which is consistent with results obtained by Dey *et al.* (2016) in the Kangra section of the Indian Himalaya where, besides the MFT, other out-of-sequence faulting such as the Jwalamukhi Thrust accommodates part of the Sub-Himalayan shortening. In our study area, one can mention either the north-propagating emerging thrust front (FBT) documented by Dasgupta *et al.* (2013) in the Brahmaputra plain, or the Main Boundary Thrust (MBT), which accommodates present-day deformation in eastern Bhutan (Marechal *et al.* 2016).

From recent studies, it is now well-established that at least two major events have occurred on the TFT in the past, the last major event having occurred about 300 yr ago (Le Roux-Mallouf *et al.* 2016; Hetényi *et al.* 2016a). Thus, a slip deficit of 3–5 m has accumulated on the TFT during this interseismic period, and could potentially be released in a large magnitude earthquake with high probability of rupture reaching the surface.

## 7 CONCLUSIONS

We have presented high-resolution near-surface geophysical imaging results based on a joint approach including electrical resistivity, seismic and gravity data to constrain the TFT geometry in south Central Bhutan.

For each data set, an inversion is performed using a stochastic approach, which combines (1) prior information on both geometric (dip angle, fault location, layer thickness) and petrophysical (resistivity, velocity and density) parameters and (2) forward models taking into account heterogeneity associated with a fault zone and superficial layers. Compared to commonly used approaches based on the search for the simplest model, the main advantages of our method include its ability to assess the fault geometry and its uncertainty on the obtained dip angle, as well as to trade-off analysis between geometric and either electrical resistivity, velocity or density properties. A joint inversion of the three data sets with physical properties linked to lithology/petrological information will be very useful, but it requires additional work that is beyond the scope of this present study.

Our results show that the upper part of the TFT is characterized by a flat and listric-ramp geometry with high variations of dip angle. This geometry clearly differs from the constant fault dip angle inferred from surface observation only. Estimating the slip rate without additional constraints from depth can therefore induce significant errors, arising both from the terrace dating process to determine the uplift rate and the projection of the fault dip angle based on surface observations.

By combining information from surface observations with our new constraints on the fault geometry, we estimate that at least 60 per cent of the Himalayan convergence is accommodated by the TFT, making this fault a zone of high seismic hazard. The hypothesis of slip partitioning cannot be totally ruled out, and other faults as the FBT emerging in the Brahmaputra plain and the MBT

can also be active. Therefore, further studies combining geomorphology and near-surface geophysics along the front, especially towards the eastern part of the Bhutan, will be useful to study potential lateral variations in the fault geometry and its implication on the present-day strain partitioning. Moreover, local variability across the TFT may be assessed by studying areas located within few hundreds of metres along strike with respect to the Sarpang study area.

## ACKNOWLEDGEMENTS

This study is funded by the French Agence Nationale de la Recherche (ANR-13-BS06-006-01). The 2015–2016 surveys benefited from the support of the French RESIF-GraviMob Research Infrastructure (http://www.resif.fr). The field campaign in southern Bhutan would not have been possible without the unfailing support of the staff and drivers of the Department of Geology and Mines, Bhutan. We sincerely extend our gratefulness to each and every one who assisted us in deployment of geophysical measurements in the field. Our thanks also go to A. Binley for providing the R2 resistivity software. We thank Mark Everett, György Hetényi and an anonymous reviewer for their helpful and constructive comments.

## REFERENCES

* Ader, T. *et al.*, 2012. Convergence rate across the Nepal Himalaya and inter-seismic coupling on the Main Himalayan Thrust: implications for seismic hazard, *J. geophys. Res.*, **117**, B04403, doi:10.1029/2011JB009071.

* Benson, A.K. & Mustoe, N.B., 1995. Analyzing shallow faulting at a site in the Wasatch fault zone, Utah, USA, by integrating seismic, gravity, magnetic, and trench data, *Eng. Geol.*, **40**(3–4), 139–156.

* Berthet, T., Ritz, J.F., Ferry, M., Pelgay, P., Cattin, R., Drukpa, D., Braucher, R. & Hetényi, G., 2014. Active tectonics of the eastern Himalaya: new constraints from the first tectonic geomorphology study in southern Bhutan, *Geology*, **42**(5), 427–430.

* Bilham, R., Larson, K. & Freymueller, J., 1997. GPS measurements of present-day convergence across the Nepal Himalaya, *Nature*, **386**(6620), 61–64.

* Binley, A., 2015. Tools and techniques: DC electrical methods, in *Treatise on Geophysics*, vol. 11, pp. 233–259, 2nd edn, ed. Schubert, G., Elsevier.

* Binley, A. & Kemna, A., 2005. Electrical methods, in *Hydrogeophysics*, pp. 129–156, eds Hubbard, S. & Rubin, Y., Springer.

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

* Cattin, R. & Avouac, J.P., 2000. Modeling mountain building and the seismic cycle in the Himalaya of Nepal, *J. geophys. Res.*, **105**(B6), 13 389–13 407.

* Cattin, R., Mazzotti, S. & Baratin, L.M., 2015. GravProcess: an easy-to-use MATLAB software to process campaign gravity data and evaluate the associated uncertainties, *Comput. Geosci.*, **81**, 20–27.

* Dahlin, T. & Zhou, B., 2004. A numerical comparison of 2D resistivity imaging with 10 electrode arrays, *Geophys. Prospect.*, **52**(5), 379–398.

* Dasgupta, S., Mazumdar, K., Moirangcha, L.H., Gupta, T.D. & Mukhopadhyay, B., 2013. Seismic landscape from Sarpang re-entrant, Bhutan Himalaya foredeep, Assam, India: constraints from geomorphology and geology, *Tectonophysics*, **592**, 130–140.

* DeCelles, P.G., Robinson, D.M. & Zandt, G., 2002. Implications of shortening in the Himalayan fold-thrust belt for uplift of the Tibetan Plateau, *Tectonics*, **21**(6), 1062, doi:10.1029/2001TC001322.

* Dey, S., Thiede, R.C., Schildgen, T.F., Wittmann, H., Bookhagen, B., Scherler, D. & Strecker, M.R., 2016. Holocene internal shortening within the northwest Sub-Himalaya: out-of-sequence faulting of the Jwalamukhi Thrust, India, *Tectonics*, **35**, 2677–2697.
Diehl, T., Singer, J., Hetényi, G., Grujic, D., Clinton, J., Giardini, D. & Kissling, E., 2017. Seismotectonics of Bhutan: Evidence for segmentation of the Eastern Himalayas and link to foreland deformation, *Earth planet. Sci. Lett.*, **471**, 2677–2697.

Dorn, C. et al., 2010. High-resolution seismic images of potentially seismic structures beneath the northwest Canterbury Plains, New Zealand, *J. geophys. Res.*, **115**, B11303, doi:10.1029/2010JB007459.

Drukpa, D., Velasco, A.A. & Doser, D.I., 2006. Seismicity in the Kingdom of Bhutan (1937–2003): evidence for crustal transcurrent deformation, *J. geophys. Res.*, **111**(6), 1–14.

Duputel, Z., Vergne, J., Rivera, L., Wittlinger, G., Farra, V. & Hetényi, G., 2016. The 2015 Gorkha earthquake: a large event illuminating the Main Himalayan Thrust fault, *Geophys. Res. Lett.*, **43**, 2517–2525.

Gautier, S., Latorre, D., Virieux, J., Deschamps, A., Skarpelos, C., Sotiriou, A., Serpersidaki, A. & Tselentis, A., 2006. A new passive tomography of the Aigion area (Gulf of Corinth, Greece) from the 2002 data set, *Pure appl. Geophys.*, **163**(2–3), 431–453.

Hammer, P. et al., 2013. Flexure of the India plate underneath the Bhutan Himalaya. *Geophys. Res. Lett.*, **40**(16), 4225–4230.

Hetényi, G., Le Roux-Mallouf, R., Berthet, T., Cattin, R., Cauzzi, C., Phuntsho, C. & Grolimund, R., 2016a, Joint approach combining damage and paleoseismology observations constrains the 1714 A.D. Bhutan earthquake at magnitude 8 ± 0.5, *Geophys. Res. Lett.*, **43**, doi:10.1002/2016GL071033.

Hetényi, G. et al., 2016b. Segmentation of the Himalayas as revealed by arc-parallel gravity anomalies, *Sci. Rep.*, **6**, 33866, doi:10.1038/srep33866.

Karastathis, V.K., Ganas, A., Makris, J., Papouilla, J., Dafnis, P., Gerolymatou, E. & Drakatos, G., 2007. The application of shallow seismic techniques in the study of active faults: the Atalanti normal fault, central Greece, *J. Appl. Geophys.*, **62**(3), 215–233.

Kumar, S., Wesnousky, S.G., Jayangondaperumal, R., Nakata, T., Kumahara, Y. & Singh, V., 2010. Paleoseismological evidence of surface faulting along the northeastern Himalayan front, India: timing, size, and spatial extent of great earthquakes, *J. geophys. Res.*, **115**, B12422, doi:10.1029/2009JB006789.

Lavè, J. & Avouac, J.P., 2000. Active folding of fluvial terraces across the Siwaliks Hills, Himalayas of central Nepal, *J. geophys. Res.*, **105**(B3), 5735–5770.

Le Roux-Mallouf, R. et al., 2015. Evidence for a wide and gently dipping Main Himalayan Thrust in western Bhutan, *Geophys. Res. Lett.*, **42**, 3257–3265.

Le Roux-Mallouf, L., Ferry, M., Ritz, J.F., Berthet, T., Cattin, R. & Drukpa, D., 2016. First paleoseismic evidence for great surface-rupturing earthquakes in the Bhutan Himalayas, *J. geophys. Res.*, **121**(10), 7271–7283.

Loke, M.H., 2015. Tutorial: 2D and 3D electrical imaging surveys, Penang, Malaysia, Universiti Sains Malaysia. Unpublished course notes, http://www.geotomosoft.com/.

Loke, M.H. & Barker, R.D., 1996. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method, *Geophys. Prospect.*, **44**, 131–152.

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

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

Molnar, P. & Tapponnier, P., 1975. Cenozoic tectonics of Asia: effects of a continental collision, *Science*, **189**(4201), 419–426.

Morey, D. & Schuster, G.T., 1999. Paleoseismicity of the Oquirrh fault, Utah, from shallow seismic tomography, *Geophys. J. Int.*, **138**, 25–35.

Mosegaard, K. & Tarantola, A., 1995. Monte Carlo sampling of solutions to inverse problems, *J. geophys. Res.*, **100**(B7), 12 431–12 447.

Nguyen, F., Garambois, S., Chardon, D., Hermitte, D., Bellier, O. & Jongmans, D., 2007. Subsurface electrical imaging of anisotropic formations affected by a slow active reverse fault, Provence, France, *J. Appl. Geophys.*, **62**(4), 338–353.

Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space, *Bull. Seismol. Soc. Am.*, **75**(4), 1135–1154.

Pandey, M.R., Tandukar, R.P., Avouac, J.P., Lave, J. & Massot, J.P., 1995. Interseismic strain accumulation on the Himalayan Crustal Ramp (Nepal), *Geophys. Res. Lett.*, **22**(7), 751–754.

Pasquet, S. et al., 2015. Detecting different water table levels in a shallow aquifer with combined P-, surface and SH-wave surveys: insights from VP/VS or Poisson's ratios, *J. Appl. Geophys.*, **113**, 38–50.

Podvin, P. & Lecomte, I., 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools, *Geophys. J. Int.*, **105**(1), 271–284.

Priolo, E. et al., 2012. The Campi Flegrei Blind Test: evaluating the imaging capability of local earthquake tomography in a volcanic area, *Int. J. Geophys.*, **2012**, 37–75.

Ramirez, A.L. et al., 2005. Stochastic inversion of electrical resistivity changes using a Markov Chain Monte Carlo approach, *J. geophys. Res.*, **110**, B02101, doi:10.1029/2004JB003449.

Schelling, D. & Arita, K., 1991. Thrust tectonics, crustal shortening, and the structure of the far-eastern Nepal Himalaya, *Tectonics*, **10**(5), 851–862.

Schuster, G.T. & Quintus-Bosz, A., 1993. Wavepath eikonal traveltime inversion: theory, *Geophysics*, **58**(9), 1314–1323.

Seeber, L. & Armbruster, J.G., 1981. Great detachment earthquakes along the Himalayan Arc and long-term forecasting, in *Earthquake Prediction*, pp. 259–277, eds Simpson, D.W. & Richards, P.G., American Geophysical Union.

Sheehan, J.R., Doll, W.E. & Mandell, W.A., 2005. An evaluation of methods and available software for seismic refraction tomography analysis, *J. Environ. Eng. Geophys.*, **10**(1), 21–34.

Singer, J., Obermann, A., Kissling, E., Fang, H., Hetényi, G. & Grujic, D., 2017. Along-strike variations in the Himalayan orogenic wedge structure in Bhutan from ambient seismic noise tomography, *Geochem. Geophys. Geosyst.*, **18**, 1483–1498.

Stephenson, W.J. & McBride, W.J., 2003. Contributions to neotectonics and seismic hazard from shallow geophysical imaging, *Tectonophysics*, **368**, 1–5.

Storz, H., Storz, W. & Jacobs, F., 2000. Electrical resistivity tomography to investigate geological structures of the Earth's upper crust, *Geophys. Prospect.*, **48**(3), 455–471.

Vernant, P. et al., 2014. Clockwise rotation of the Brahmaputra Valley relative to India: tectonic convergence in the eastern Himalaya, Naga Hills and Shillong Plateau, *J. geophys. Res.*, **119**, 6558–6571.

Villani, F., Tulliani, V., Sapia, V., Fierro, E., Civico, R. & Pantosti, D., 2015. 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**(3), 1482–1494.

Wise, D.J., Cassidy, J. & Locke, C.A., 2003. Geophysical imaging of the Quaternary Wairoa North Fault, New Zealand: a case study, *J. Appl. Geophys.*, **53**(1), 1–16.

Won, I.J. & Bevis, M., 1987. Computing the gravitational and magnetic anomalies due to a polygon: Algorithms and Fortran subroutines, *Geo-physics*, **52**(2), 232–238.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Bivariate frequency histograms between dip angle in degrees and the other parameters of the ERT model obtained from Wenner–Schlumberger array and electrode spacing of 1 m. Colour scale gives the probability of occurrence of parameter pairs within our a priori range of parameters.

**Figure S2.** Same as Fig. S1, except for electrode spacing of 2.5 m.

**Figure S3.** Same as Fig. S1, except for electrode spacing of 5 m.

**Figure S4.** Relationship between the obtained dip angle in degrees and the other parameters of the velocity model.

**Figure S5.** Relationship between the obtained dip angle in degrees and the other parameters of the density model. Note that density
means density contrast with respect to density of the north deep layer (NL).

**Figure S6.** Res2dinv-generated inverse model resistivity section for Wenner–Schlumberger array with 5 m spacing (bottom). The observed (top) and calculated (middle) apparent resistivity pseudo-sections are also shown.

**Figure S7.** Top: tomographic image showing the velocity variations of both sides of the TFT using a threshold value of 100 for the ray coverage. After 20 iterations, we obtained traveltime residual of 1.41 ms corresponding to an rms value of 3.1 per cent. Bottom: ray coverage illustrating the resolved area.

**Table S1.** Inversion results. The range column gives the prior pa- rameter range used in the inversion. The obtained parameter values are associated with the highest relative frequency model. The asso- ciated uncertainties are indicated in brackets. Uncertainties on each parameter are given by the full width at half maximum. The symbol '-' means no constraint has been obtained.

Please note: Oxford University Press is not responsible for the con- tent or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be di- rected to the corresponding author for the paper.