# Shallow subsurface morphotectonics of the NW Sumatra subduction system using an integrated seismic imaging technique

Dibakar Ghosal, S. C. Singh, J. Martin

► To cite this version:

Dibakar Ghosal, S. C. Singh, J. Martin. Shallow subsurface morphotectonics of the NW Sumatra subduction system using an integrated seismic imaging technique. Geophysical Journal International, 2014, 198, pp.1818-1831. 10.1093/gji/ggu182. insu-03581103

HAL Id: insu-03581103

https://insu.hal.science/insu-03581103

Submitted on 19 Feb 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Shallow subsurface morphotectonics of the NW Sumatra subduction system using an integrated seismic imaging technique

Dibakar Ghosal,* S. C. Singh and J. Martin

Laboratoire de Géosciences Marines, Institut de Physique du Globe de Paris, 1 rue Jussieu, Cedex-05, F-75238 Paris, France.

E-mail: dibakarghosal@iitbbs.ac.in

Accepted 2014 May 15. Received 2014 May 15; in original form 2013 August 15

## SUMMARY

In subduction zones, shallow subsurface structures are the manifestation of the plate interactions at depth. However, significant water depths, rough bathymetry and presence of heavily deformed accretionary wedge materials hamper imaging of the near-surface features to a great extent using conventional imaging techniques. In this study, we show results using an integrated processing technique to a multichannel seismic data set acquired in 2006 from the northwestern offshore Sumatra. We start with first downward continuing the 12-km-long surface streamer data to the seafloor, followed by a high-resolution traveltime tomography of refracted phases to determine a detailed velocity-depth model of subsurface, which in turns, is used for pre-stack depth migration in order to delineate the shallow subsurface structures beneath the trench, subduction front and outer accretionary wedge. Our velocity-depth model and the depth migrated image depict variation of sediment properties across the front and structures of uppermost sedimentary sequence with an unprecedented high resolution providing the precise location of the frontal and conjugate thrusts, highly folded sedimentary sequences, which in turns describe their relationship with the top of the subducting plate and factors that control rupture propagation to the trench. Furthermore, we estimate the porosity distribution across the front, where we find a 12 and 18 per cent decrease in porosity beneath the deformation front and the inner accretionary plateau at 500 m below the seafloor, respectively, which we interpret to be due to the compaction. A significant decrease in porosity at the plate interface below 5–6-km thick sediments near the deformation front would increase the coupling, leading to the rupture propagation up to the trench, uplifting 4.5 km water and producing large tsunami.

**Key words:** Controlled source seismology; Seismic tomography; Subduction zone processes; Continental margins: convergent; Folds and folding.

## 1 INTRODUCTION

Accretionary wedges associated with convergent margins are developed by the accretion of scraped off sediments from the top of a subducting oceanic plate against a continental backstop or island arcs (Byrne et al. 1988; Saffer & Bekins 2002). Accretion of such water-rich unconsolidated or semi-consolidated materials cause the outer wedge (accretionary prism) to slip aseismically (Byrne et al. 1988), whereas the inner wedge (forearc), which overlies the seismogenic zone, is usually comprised of old consolidated materials and slips seismically causing great megathrust earthquakes in the subduction zones. Wedge dynamics is a complicated process and has been investigated vigorously since the early 1980s (Dahlen et al. 1984; Dahlen 1990; Wang & Hu 2006). The morphology and the

characteristics of the wedges forming materials are well constrained at many convergent margins (e.g. Nankai, Barbados, Cascadia) from high-resolution bathymetry, seismic and borehole data sets (Von Heune & Scholl 1991; Bangs et al. 1996; Calvert 2004; Calvert et al. 2006). However, it is poorly constrained along the western offshore Sumatra, which has produced three major earthquakes of $M_w > 8.5$ in the last decade. Our present study is focused on the NW Sumatra margin that contains a 160–180-km wide accretionary wedge (Singh et al. 2008) with highly deformed sedimentary landforms associated with numerous folds, bivergent faults and fractures.

In order to image the surface morphology and sedimentary structure of the North Sumatra subduction zone, several bathymetry, seismic reflection and refraction surveys were carried out in the recent years in the epicentral region of the 2004 December 26 earthquake (Henstock et al. 2006; Franke et al. 2008; Graindorge et al. 2008; Singh et al. 2008). Seismic images presented so far have been either in the time domain or converted into depth using a large-scale velocity structure determined from ocean bottom

<footnote>* Now at: Indian Institute of Technology Bhubaneswar, Toshali Bhawan, Satyanagar, Orissa, India.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Tectonic map of Sumatra subducting system. GSF, Great Sumatra Fault, WAF, West Andaman Fault, MT, Mentawai Fault, WSC, Wharton Spreading Centre. Different major earthquake events are shown by red stars. D-AA' and S-22 are profiles used in Dean *et al.* (2010) and Gulick *et al.* (2011). Yellow rectangle marks region shown in fig. 1b. F6, F5—Fractures 6 and 5, respectively. (b) Detailed study area. Bathymetry data were compiled from the French surveys (Singh *et al.* 2008) with GEBCO data set in background. Different striking features are marked. Dashed black line indicates the extrapolation of fracture zone F6 (Singh *et al.* 2011a). North is marked by the arrow.</figcaption></figure>

seismometer (OBS) data (Franke et al. 2008; Singh et al. 2012), but to study the deformation front a detailed high-resolution depth image of the subduction zone is required.

Here, we present results from a high-resolution seismic reflection, refraction and bathymetry survey which were carried out in 2006 at the 250 km NW from the epicentre of 2004 December 26 earthquake (Fig. 1). In order to determine the shallow subsurface velocity, we apply a high-resolution traveltime tomography (TTT) on downward continued shot and receiver gathers, and then we use this velocity-depth model for depth migration of seismic reflection data. Furthermore, we use this high-resolution velocity-depth model to estimate the porosity distribution along the profile. All these results are used to shed light upon the characteristics of sedimentary structures across the deformation front and their influence on large earthquakes at the offshore northern Sumatra.

## 2 STUDY AREA

Northwest offshore Sumatra contains a sediment-filled trench, a striking deformation front and a wide conspicuous accretionary wedge (Fig. 1b), and is surrounded by the Nicobar Island in the north, fossil Wharton spreading centre in the south, the Aceh forearc basin in the northeast and the prolonged Ninety-East ridge in the west. The trench in this area consists of 2–5-km thick sedimentary deposits overlying a 58–60-Ma-old oceanic crust (Liu et al. 1983) formed at the Wharton spreading centre system (Figs 1a and b), and several fracture zones (Singh et al. 2011a) that have been reactivated as sinistral strike-slip faults (Deplus et al. 1998; Carton et al. 2014)

due to the diffuse deformation in the Central Indian Ocean (Gordon et al. 1998). The terrigenous sediments, originated either by the denudation of the vast Himalayan mountain belt or by the erosion of the backarc structural highs, started depositing in Late Oligocene (Matson & Moore 1992) which are, subsequently transported by fluvial mechanisms and accreted through subduction processes against the continental hard and crystalline backstop, leading to the formation of the wide accretionary wedge. In the accretionary wedge, numerous NW–SE trending fold-thrust belts and highly deformed structures have developed due to the trench normal thrust component associated with the oblique subduction of the downgoing plate. The rate of subduction in this area lies between 52 and 47 mm yr⁻¹, with increasing obliquity towards west (Fig. 1a). The present study area lie in the 2004 December 26 earthquake (M<sub>w</sub> = 9.3) rupture zone.

## 3 DATA ACQUISITION

During 2006 July–August, coincident deep seismic refraction and reflection profiles were shot using the French R/V Marion Dufresne and the WesternGeco M/V Geco Searcher vessels carrying 8260 and 10 170 cubic inch airgun array sources, respectively. The WG2 profile (Fig. 1b) is oriented 20° anticlockwise from the trench normal on which 56 OBS spaced at 8.1 km were deployed and shots were fired at 150-m intervals (Chauhan 2010). A 12-km-long Qmarine streamer containing 3832 individual hydrophones spaced at 3.125 m was deployed at 15-m water depth, which was then grouped at 12.5 m after a digital filter (Martin et al. 2000), record deep seismic reflection energy at 50-m shot interval (Singh et al. 2012).
Another 5.5-km-long Q-marine streamer was towed at 7.5-m water depth to acquire high-resolution reflection data for better imaging the near-surface features. The WG2 profile covers the trench de-posits, accretionary wedge, forearc basin, Sumatra platform and arc-related basins. The forearc and the arc-related basins are in-vestigated by Singh *et al.* (2012) and Ghosal *et al.* (2012); in this study we focus on the 200-km segment that covers a part of the trench—and the outer accretionary wedge (Fig. 1b).

# 4 METHODOLOGIES AND RESULTS

## 4.1 Bathymetry around WG2 profile

The bathymetry data show that water depth varies drastically from the trench at 4.5–0.3 km on the accretionary wedge. The accre-tionary prism contains four ridge-like structures marked as R1, R2, R3 and R4 (Fig. 2), which are bounded by thrusts and separated from each other by piggyback basins. The seafloor is nearly flat (~4.5 km) in the trench (Fig. 2a) up to the subduction front, the toe of the ac-cretionary prism, whereas the water depth decreases updip along the frontal slope (5°–6°) and reduces to 337.5 m on top of the frontal ridge R4 (Fig. 2b). The water depth suddenly increases to 1600 m on the northeast side of the frontal ridge and gradually decreases along the gently varying accretionary plateau (Fig. 2b), which consists of NW–SE oriented fold-thrust belts containing small ridge-like structures. The accretionary plateau is subdivided into two parts: outer plateau and inner plateau depending upon its surface slope that ranges between 1° and 3°. The outer plateau is extended up to nearly 55 km and is separated from the inner plateau by a piggyback basin.

## 4.2 High-resolution P-wave velocity estimation

Despite having a long streamer (12 km), conventional techniques of high-resolution velocity estimation are not well suited to this part of the study area due to several interrelated reasons: (1) strong dips on the seafloor, folding and faulting of the sediments and rough

seafloor topography do not permit accurate velocity estimation using the normal moveout velocity method, (2) the highly deformed and compacted sediments do not produce coherent reflections, par-ticularly on the plateau and (3) the absence of refraction arrivals in deep water (Fig. 3a). The 8.1-km OBS spacing was too large to determine detailed near-surface velocity. In order to enhance the refraction arrivals and bring them before the seafloor reflections (Arnulf *et al.* 2011), we have downward continued both shots and receivers onto the seafloor and have carried out first arrival TTT.

### 4.2.1 Downward continuation (DC)

The DC is an efficient technique (Berryhill 1979; Arnulf et al. 2011) to extrapolate the recorded wavefield to an arbitrary surface. Apart from bringing the refracted arrival as first arrivals, the DC method collapses seafloor diffractions thus improving the imaging condition (Arnulf et al. 2011). Practically, DC can be carried out in different ways: (1) in the frequency domain—using pre-stack phase-shift method (Gazdag 1978), (2) in the time domain (Berryhill 1979)—applying Kirchhoff’s integral method and (3) solving the wave equation using finite difference method (Claerbout 1976).

Here, we have used the Kirchhoff’s integral extrapolation tech-nique in the time domain which is formulated by convolving the recorded wavefield at the water surface with the time lagged delta functions and summing these convolved results all over the receiver locations for all shot gathers (Harding et al. 2007; Fig. S1). The time lags are estimated by dividing the distance between the sur-face locations of the shot gathers and their extrapolated positions on seafloor by a constant water velocity. Since no salinity log along the profile was available we used a constant water velocity of 1480 m s⁻¹. The DC process is carried out in three steps: (1) extrapolation of the receiver gathers vertically down to the seafloor, (2) sorting of the extrapolated shot gathers into common receiver location gathers (CRL) and (3) DC of the shot positions to the seafloor and sorting back into shot gathers. The aliasing effect is reduced by interpo-lation of shots to 12.5 m and applying a bandpass frequency filter of 5–18 Hz after the step 2. Fig. 3(c) shows the results of DC of

<figure><img src="image_3.png" /><figcaption>**Figure 2.** 3-D perspective view of the seafloor bathymetry. (a) Frontal slope: marks different ridges (R1, R2) and WG2 profile; (b) Accretionary wedge: white lines indicate the NW–SE lineations caused by the trench normal component of oblique subduction; AP, Accretionary plateau.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Surface 12-km-long streamer filtered shot record and associated ray diagram (a, b) for surface seismic and the same shot record after downward continuation (DC) to place the source and receiver positions on the seafloor (c, d). The shot (Sh) and streamer positions (R1–R956) are marked by red star and the yellow line respectively.</figcaption></figure>

data near the deformation front. On the surface seismic data, the reflection part of the data is contaminated by seafloor scattering effect (Fig. 3a) and seismic refraction can only be observed at offsets from 9 to 12 km (Figs 3a and b). In the downward continued data, the refraction arrivals can be observed from 0 to 9 km offset. Figs 3(b) and (d) show the ray diagram for turning rays for surface seismic and downward continued shot gathers, clearly indicating that the downward continued data have much higher ray coverage, particularly for the shallow part.

### 4.2.2 Traveltime tomography

The TTT was performed by inverting the first arrivals (refracted) of every fourth downward continued shot gathers to determine a velocity–depth model using the algorithm developed by Van Avendonk et al. (2004). Traveltimes are picked using a semi-automatic technique, which interpolates traveltimes between two picked traveltimes and reduce the picking time. Consistency in identifying the first arrivals between two consecutive shot gathers are maintained using the preceding traveltime curves as a guide function for the following picks. In total, 600 000 of traveltimes were picked. The uncertainty assigned to each of the picks are based on several factors such as: (1) the shot and receiver positions on the water surface lead to a traveltime uncertainty of 2 ms, (2) uncertainty associated with the bathymetry is set to ~10 ms, (3) the uncertainty linked to the shot and receiver positions on the seafloor after DC may lead to an approximate traveltime uncertainty of 8 ms and (4) uncertainty associated with the traveltime picking is estimated to be 8 ms. Considering these uncertainties as mutually independent, we determine the overall uncertainty as 15 ms.

The traveltimes were computed using the shortest path method (Moser 1991) in a model, which was discretized into regular 50-m square grids. The starting model was based on the pre-existing 2-D velocity–depth model (Fig. 4a) obtained through the tomography results of the OBS experiment conducted coincidentally along the WG2 profile (Chauhan et al. 2009; Singh et al. 2012). The ray diagram in the starting model shows uniform ray coverage down to 3 km below the seafloor (bsf) at this location. The misfit between the observed traveltime and computed traveltime in the starting model varies from ~50 ms at the farthest offset to ~150 ms at near offsets, suggesting that there is a significant mismatch between the starting model and real model in the uppermost part of the model.

The misfit between the observed and calculated traveltime is minimized in a least-square sense. The regularization of the inversion is controlled by a roughness matrix, which contains first and second-order derivatives. Assuming the velocity is rougher in the vertical direction, the horizontal derivatives were given four times higher weight than that of the vertical one during the inversion, whereas the model artefacts associated with the inversion are suppressed by proper choice of the damping constant (Van Avendonk et al. 2004). The overall misfit is measured by $\chi^2$ which is the sum of the squared individual traveltime misfits divided by the overall pick uncertainty. The misfit is gradually reduced through a series of linearized inversions until the desired fit is achieved, typically $\chi^2 = 1$, at which point the average misfit matches the pick uncertainty. At each iteration, velocity nodes are updated throughout the ray covered portion of the model.

#### 4.2.2.1 Inversion results
The best-fitting P-wave velocity–depth model is shown in Fig. 4(b). After 18th iteration, the root mean
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Inversion results. (a) Starting model adapted from the OBS tomography result (Chauhan 2010) (b) Best-fitting model after 18th iteration. Dotted black lines indicate velocity contours. (c) Interpreted image describes different layers associated with the subsurface. Three major layers are identified. L1: Unconsolidated layer, L2: Semi-consolidated to consolidated layer and L3: Highly compacted and lithified basement type deposits. CT, Conjugate thrust, FT, Frontal thrust and different ridges are marked as R1, R2, R3 and R4.</figcaption></figure>

square misfit for the best-fitting model reached to 15 ms (Fig. 5), which is the same to the pre-defined uncertainty in the data picking. The ray coverage can be visualized using the derivative weight sum (DWS) plot (Fig. 6), which is derived by summing the columns of the Fréchet derivative matrix (Scales 1987). Although this plot does not provide actual information on the resolution, it does give information about the ray coverage or the sensitivity of a model parameter to the data (Van Avendonk et al. 2001). Traced rays sample the uppermost sedimentary layers down to a depth of 2.5–3 km bsf beneath the trench, whereas the ray coverage is deeper (3.5–4 km bsf on an average) in the accretionary wedge (Fig. 6).

The velocity model (Fig. 4b) within the trench shows that the ve-
locity at the topmost semi-horizontal sedimentary layers lie between
1.6 and 3.0 km s⁻¹. Further northeast within the accretionary prism,
the velocity contours run updip following the shape of the anticlinal
ridges. The 1-D vertical velocity profiles (Fig. 7), extracted from the
2-D velocity–depth model (Fig. 4b), show two distinct slope breaks
along the profiles; a layer containing high-velocity gradient over-
lies a layer with low-velocity gradient followed by a layer of high-
velocity gradient. Based on the changes in slopes of the 1-D profiles,
the overlying sediments are subdivided into three layers with average

velocities: (1) L1 (1.5–2.7 km s⁻¹), (2) L2 (2.7–3.7 km s⁻¹) and (3) L3 (>3.7 km s⁻¹) (Fig. 4c). The sediment deposits (~0.5–1.5 km) in layer L1, corresponds to the first high-velocity gradient, suggesting the presence of the uncompacted sediments (Fig. 4c), whereas the underlying L2 layer with low-velocity gradient implies the existence of a semi-compacted to compacted sediments (Fig. 4c). The L3 layer is identified as the second high-velocity gradient at the bottom of 1-D profiles, representing the transition from the compacted high-velocity layer deposits to the underlying highly compacted or crystalline sediments (Fig. 4c).

We suggest that increase in velocity is related to compaction and release of fluids from the sediments and thus compaction increases from L1 to L3. Beneath R4, the 2.5 km s⁻¹ velocity contour touches the seafloor, suggesting the presence of severe folding and subsequent erosion leading to the increase in velocity beneath the ridge. Beneath the accretionary plateau, the velocity contours extend northeastwards following the bathymetry and show severe undulations.

From the velocity contours, one can clearly observe that the
layer L1 is thick in the trench fill area and thins up northeastwards
along the frontal slope, whereas layers L2 and L3 are observed
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Traveltime misfit plots: (a) Starting misfit (red), final misfit (blue); (b) Histogram of the data errors for all shots and receivers. A total of 600 000 traveltimes were picked. The shape of the distribution is nearly the Gaussian. The width of the bars is 15 ms.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** The derivative weight sum (DWS) for the best-fitting model. The DWS is calculated adding the elements of each column of the Fréchet derivative matrix. It is a dimensionless measure of the sensitivity of the data to a model parameter (Van Avendonk et al. 1998). The DWS is plotted between 0.1 and 10. Ray coverage is maximum below the frontal slope and accretionary wedge down to 1.5–2 km bsf.</figcaption></figure>

predominantly beneath the outer accretionary wedge. The base of the layer L3 has velocity >4.5 km s⁻¹ and is observed at a depth of 3.5 km bsf suggesting the existence of a highly compacted to crystalline-like sediments at depth.

## 4.3 Pre-stack depth migration results

Both the long and short streamer data sets are used for the pre-stack depth migration. The interval velocity was obtained by combining the above high-resolution velocity with the low-resolution OBS tomography (Chauhan 2010) velocity using a linear interpolation down to 12 km depth in order to image the top of the subducting oceanic crust below the accretionary prism and its relation with the overlying shallow structures (Fig. S3).

Some pre-processing including swell noise attenuation and multiple removal were performed before starting the migration. The pre-processing of the long streamer data was carried out at WesternGeco, in Jakarta (Singh et al. 2012) using 8-ms sampling interval

thus making the data sets very useful in delineating the deeper structures. For the short streamer, we kept the sampling interval to 2 ms in order to obtain high-resolution image. The swell noise was attenuated using a filter in frequency–wavenumber (FK) domain (Sheriff & Geldart 1982). Subsequently, we suppressed the multiples using a radon multiple removal technique (Sheriff & Geldart 1982), which efficiently eliminates the multiples from the far offset sections. The residual multiples from the near offset are effectively suppressed by applying an inner trace mute (Fig. S4). The side scattering effects were reduced using a dip filtering technique. During the migration, the aperture length is set to 3.5 km and a mild triangular antialiased filter is applied.

The uppermost part of the trench deposits is composed of north-eastward gently dipping thick turbiditic sediments, whose thickness increases from 2.5 km in the southwest to 5.5 km near the subduction front (Figs 8 and 9). Sedimentary layers are transected by numerous vertical to subvertical normal faults with dips between 70° and 80°, mainly dipping towards NE. Most of these normal faults stop within the sedimentary layers and do not arrive to the seafloor,
<figure><img src="image_9.png" /><figcaption>**Figure 7.** 1-D vertical velocity sections extracted from Fig. 4(b) at different locations along the profile. (Red) Seafloor corrected (SFC) 1-D velocity for OBS tomography results and (Blue) Seafloor corrected (SFC) 1-D velocity for high-resolution tomography results. Different layers L1, L2 and L3 are marked by the black arrows.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Pre-stack depth migrated image of the long streamer data using a combined interval velocities obtained from the high-resolution tomography and the OBS tomography. (a) Uninterpreted and (b) Interpreted seismic images. Blue: Basement, Purple: High-amplitude Reflection Packets (HARPs) and Red: Different faults. CT, Conjugate thrust, FT, Frontal thrust, T1, Landward-dipping thrusts, various ridges on the frontal slope are marked by R1, R2, R3 and R4. S1 and S2 are slopes of the bathymetry on the accretionary plateau. (c) Combined velocity model is superimposed on the interpreted depth section and velocity contours are marked by thin black lines. Regions enclosed by the black and blue rectangles are described in Figs 9 and 10, respectively; vertical dashed lines are positions of 1-D velocity profiles described in Fig. 11. The 40 and 4 per cent porosity contours are marked by the white lines.</figcaption></figure>

suggesting that they are not active, but few of them, particularly between 0 and 20 km distance range, reach to the seafloor, suggesting that they are active, and likely to be formed either by plate bending effect or dewatering (Fig. 9). Some weakly reflective, localized structures are commonly found within the sedimentary layers, which might have resulted from gas accumulation due to the lithification and subsequent degassing from the deeper layer or could be buried palaeochannels, mostly originated from the adjacent frontal deposits. These types of channel deposits are commonly observed

in other parts of offshore northern Sumatra (Mosher *et al.* 2008) and in some other subduction environments as well (Huyghe *et al.* 2004; Goldfinger & McNeill 2006) representing mass wasting associated with the frontal slope. Near the base of the sedimentary layers, a zone of high-amplitude reflection packets (HARPs) with a thickness of 300 m is observed with irregular shape which indicates substantial contrasts in impedance (Fig. 9). A similar layer is observed further south (Fig. 1) and is interpreted as a pre-decollement layer—a seaward extension of rupture propagating surface for megathrust (Dean
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Detailed image (scale 1:1) of the trench deposits using the pre-stack depth migration of the long streamer data sets. Red lines are faults in the trench deposits; dotted circle indicates evidence of palaeochannel; blue and black arrows indicate the HARP and top of the subducting plate, respectively; the polarity of the high-amplitude reflector is compared with the seafloor and top of the oceanic crust from the region enclosed in the cyan rectangle. σ max: maximum principal stress.</figcaption></figure>

et al. 2010). Since this layer has a ~300-m thickness and irregular in shape, we interpret this as a lithologic boundary that might have formed due to deposition of pelagic sediments (Singh et al. 2008, 2011b) or thin beds of sand silt alterations, which are often observed in alluvium fans (Curray et al. 2003). A low-amplitude zone with average thickness of 1 km is sandwiched between the HARP and the underlying gently dipping (ranges from 2° to 3°) subducting plate, which is well imaged between 8 and 10.5 km depth within the trench. The top of the downgoing slab (basement) is highly irregular and displays some elevated ridge-like structures at 15 and 20 km distance ranges SW of the trench. Profile WG2 crosses a fracture zone (F6, in Singh et al. 2011a) and these basement features are presumably due to the fracture zone topography.

The 1-D velocity–depth (Fig. 11) sections extracted from the combined velocity–depth (Fig. 8c) model illustrate the velocity variation between the HARPs and the top of the oceanic crust (TOP) inside the trench deposits. The velocity–depth profile suggests that the velocity estimated between HARP and TOP varies from ~5.1 to 5.5 km s⁻¹ close to the front, whereas at 27 km offset (middle of the trench) it ranges between ~4.8 and 5.5 km s⁻¹. The velocity at the base of the sedimentary layer is very high suggesting that the sediments are highly compacted causing strong locking with the subducting oceanic slab and develops a suitable regime for the rupture to propagate from the seismogenic zone to the front and further seawards. We do not observe any low velocity below HARP as indicated by Dean et al. (2010), but we do observe low-amplitude reflectivity below the HARP.

The gently dipping layered sedimentary sequences of the trench deposits start to get deformed at the subduction front, which is located at the toe of the accretionary wedge, producing severe folding and faulting inside the prism (Fig. 10). The frontal anticlinal ridges (R1 and R2) are mostly bounded by bivergent fault systems. Due to significant amount of compaction during the ongoing accretion related to the subduction process, the degree of folding of the sedimentary layers inside these ridges increases updip along the frontal

slope and attains a steepness of 40° in the ridge R4. The dip of these structures are estimated accurately, and while the frontal slope is estimated as 6°, the surface slopes of ridges R1, R2, R3 and R4 are 5°, 9°, 10° and 10°, respectively (Fig. 10). The prism taper is estimated as 10° combining the frontal slope (6°) with the dip of the basement (4°) below the frontal slope and suggests that the prism taper is high, which in turn implies that the basal friction is strong (Dahlen et al. 1984). The bivergent faults surrounding the ridges show relatively steep slope. The frontal ridge (R1) is 11-km wide and is bounded by steeply dipping bivergent thrusts (dip ~ 43°) that seem to join within the trench fill sediments. The cumulative slip (Hill et al. 2012) on the landward-dipping branch is about 100 m whereas that on the seaward-dipping branch is nearly five times of that (550 m), suggesting that the most of the slip is taken up by seaward-dipping branch. The landward-dipping branch seems to continue down to the basement and joins the decollement surface at a triple junction. A shallow landward-dipping fault lies beneath R1, just above the basement. A large slip seems to have occur along the CT3 (4 km), suggesting that the triple junction might be the seaward tip of the megathrust, producing a pop-up structure near the subduction front, responsible for a large uplift and large tsunami. The top of the subducting oceanic slab is poorly imaged throughout the accretionary wedge, but a small segment of it is observed at the toe of the prism (Fig. 10). The ridges R3 and R4 are bounded by landward-dipping thrusts. All the thrusts SE of T1 are seaward-dipping whereas the FT-CT3 system shows different folding and faulting pattern. The horizontal layers beneath thrust T1 and CT3 and the folding and faulting between FT and CT3 suggest T1 was active for a long time, and then jumped to the present position in the last few hundred thousand years (considering total slip due to thrust (T1) as 4.5 km and the rate of convergence as 45 mm yr⁻¹).

The NE end of the frontal slope (at 100-km distance) is demarcated by a sudden drop in slope on the accretionary plateau, which is ~100-km wide and comprised of two major structural elements: outer plateau and inner plateau with the average slopes of 1° (S1)
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Detailed image (scale 1:1) below the frontal slope using the pre-stack depth migration of the short streamer data sets. (a) Uninterpreted (b) Interpreted. Red: various faults associated with the ridges. Frontal thrust (FT), the conjugate thrusts (CT) and landward-dipping thrusts (T1) are marked by the red lines. Sedimentary beds are marked by the green lines. σ max: maximum principal stress.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** 1-D NMO (red) and combined (blue) velocity-depth models are shown at 27 and 55 km offsets in the trench deposits. HAR, high-amplitude reflector and TOP, top of the oceanic plate. Combined velocity model shows similarity with the NMO at the uppermost section and becomes higher at depth.</figcaption></figure>

and 3° (S2) (Fig. 8). The outer plateau (S1) is comprised of five symmetric ridges, each bounded by thrust faults and small piggyback basins and a wide basin that contains thick (3.5 km) unconsolidated sediments (similar to the frontal sediments) that are deformed internally and starkly differ from the neighbouring compacted sediments. These sediments must have resulted from the jump in the position of the main thrust as mentioned above. A 20-km wide asymmetric basin, which lies at 1700 m water depth, separates the two plateaus and is filled with subhorizontal sedimentary deposits. The inner plateau forms the forearc high, and is comprised of poorly reflective material. Chauhan (2010) suggests that this forearc high is bounded in the northeast by backthrusts separating the forearc high with the Aceh forearc basin.

## 4.4 Porosity-depth estimation

Characterization of the porosity-depth variation underneath the trench deposits and accretionary wedge is computed using the velocity derived from the high-resolution TTT. Since no borehole data sets are available along the offshore northern Sumatra margin, we computed the porosity using an effective medium theory (EMT; O'Connell & Budiansky 1974; Berryman 1980; Norris 1985; Jakobsen et al. 2000). Subsequently, the results are compared with other subduction zones such as Nankai and Barbados, where the borehole data sets are available down to a depth <1 km bsf (Gamage et al. 2011).

The EMT, which is more promising than any simple averaging methods (Hill 1952; Hashin & Shtrikman 1963), is used to relate the crack and pore properties with the seismic properties following two different approaches: differential effective medium (DEM) theory (Norris 1985; Berge et al. 1992; Berryman 2002) and self-consistent approximations (SCAs; O'Connell & Budiansky 1974). The DEM is useful to explain a composite that contains isolated cracks comparable to the deposits above the decollement surface (Tsuji et al. 2008) in the subduction zones and for igneous crust at the mid-oceanic ridge systems (Christensen 1984; Taylor & Singh 2002; Seher et al. 2010). The crack geometry is defined by the aspect ratio (r), which is in turn determined by the ratio of the polar and equatorial radius of inclusions. The shape and orientation of these inclusions are very important in determining the nature of the predicted velocities. The predicted velocity would be isotropic if the inclusions are elliptical or spherical and oriented randomly, whereas it would be anisotropic for the aligned elliptical intrusions (Taylor & Singh 2002). Moreover, the crack orientations are responsible for producing fast and slow velocities for a composite at a particular aspect ratio.

We assume a two-phase EMT that requires a background matrix and a fluid phase. The background matrix is composed of four other constituents: sand, silt, calcite and clay. The fluid phase is considered
<table><caption>**Table 1.** Bulk modulus and shear modulus of different composites are listed.</caption><thead><tr><th>Composites</th><th>Bulk modulus (k) (GPa)</th><th>Shear modulus (μ) (GPa)</th></tr></thead><tbody><tr><td>Sand</td><td>37.79</td><td>44.07</td></tr><tr><td>Silt</td><td>37.645</td><td>29.535</td></tr><tr><td>Clay</td><td>21.2</td><td>6.66</td></tr><tr><td>Calcite</td><td>76.8</td><td>32.0</td></tr><tr><td>Seawater</td><td>2.25</td><td>0.0</td></tr></tbody></table>

as sea water (Table 1). We have implemented the Voigt–Reuss–Hill (VRH) averaging method (Hill 1952) to calculate the effective bulk modulus and shear modulus for the background matrix. The VRH averaging technique calculates an arithmetic mean of the Voigt and Reuss estimates for isotropic aggregates based on their elastic properties of the matrix materials. These effective elastic moduli are used as input for the DEM approach for different aspect ratios ($r = 1–100$) (Jakobsen et al. 2000). The P-wave velocity of the pure matrix is 4.8 km s$^{-1}$. Throughout the calculations, the effects of different intrinsic parameters (e.g. temperature and permeability) for the medium are ignored.

We start the DEM calculation assuming the existence of the solid background matrix with 0 per cent porosity and eventually we increase the porosity until a pre-defined value. Experimental data suggest that the porosity of the sedimentary deposits lies in a range of 50–70 per cent for a few hundred metres below the seafloor (Gamage et al. 2011) and becomes <30 per cent just above the decollement surface (Tsuji et al. 2008). Therefore, we have set the maximum porosity value as 70 per cent. We compare the DEM results for different aspect ratios ($r = 1, 10$ and $100$) with the borehole data sets observed at two different subduction zones: Nankai Trough and Barbados. Since subduction parameters for the Nankai Trough (rate of subduction is 4 cm yr$^{-1}$ and sedimentary thickness on accretionary prism is 1.0–2.0 km) is more similar to the NW Sumatra subduction zone, we have fixed the partial volumes of different matrix building constituents (Gamage et al. 2011) close to the parameters observed at Nankai Trough.

We first estimate 1-D porosity for the trench deposits at 20 km SW of the subduction front for different aspect ratios and compare these with that observed from Barbados and Nankai Trough (Fig. 12c). The porosity–depth results display least fit with the borehole data sets for spherical inclusions ($r = 1$), whereas randomly oriented elliptical inclusions ($r = 100$) better explain the observed data (Fig. 12c). Given the fit between the borehole data and our model, we use $r = 100$ and randomly oriented inclusions for further discussion. The velocity–depth model at 75-km distance indicates an increase in velocity by 55 per cent at 2 km depth bsf, whereas the porosity–depth model shows a decrease in porosity by 60 per cent. It is interesting to note that the porosity rapidly decreases from 70 per cent at the seafloor to 12 per cent at 1.7 km bsf, and further decreases from 12 per cent to 3 per cent between 1.7 and 3.3 km depth, suggesting that the most of the fluids are expulsed in the upper 1.7 km bsf, and the sediments below this depth are presumably well compacted. Fig. 8 shows that the porosity contours dip opposite to the dip of the trench sediments indicating that the porosity reduces laterally landwards within the trench, suggesting that the compaction associated fluid expulsion initiates 40 km seawards from the deformation front, implying the decollement to extend 40 km seawards. Considering the average porosity of the trench deposits as 17 per cent at 1 km depth bsf at 45 km distance, we have estimated the average porosity at 1 km bsf on the subduction front

<figure><img src="image_14.png" /><figcaption>**Figure 12.** (a) Porosity–velocity relationship established using the differential effective medium (DEM) theory for disordered crack inclusions for different aspect ratios ($r = 1, 10, 100$). (b) 1-D velocity–depth model extracted from the high-resolution tomography results Fig. 4(b) at the 25 km NE from the trench. (c) Porosity–depth relationship for different aspect ratios of randomly oriented cracks determined using the velocity–depth result from the tomography. Porosity–depth results for the real data sets from the Nankai (black circles) and Barbados (blue stars) subduction regions.</figcaption></figure>

(80 km distance), on the outer accretionary plateau (135 km distance), and on the inner accretionary plateau (180 km distance) and found them to be 16, 13 and 11 per cent, respectively. The details of the porosity distribution at 0.5, 1 and 1.5 km bsf are explained in Table 2. These results suggest that porosity decreases arcwards from the trench indicating the lateral compression associated with the subduction process. It also decreases with depth as well, requiring the increase of overburden stress associated with thickening of sedimentary layers.

## 5 DISCUSSION

The 2.0–5.0 km thick, highly reflective, turbiditic deposits, which are primarily derived from the Bengal and Nicobar fans, have been deposited since early to mid-Eocene time (Curray & Moore 1971)
<table><caption>**Table 2.** Porosity decreases with depth and landwards due to lithification and lateral compaction, respectively. Negative signs indicate the percentage of porosity decrease with respect to trench deposits at different depths.</caption><thead><tr><th>Depth bsf (km)</th><th>Porosity at trench (per cent)</th><th>Frontal slope (per cent)</th><th>Outer accretionary plateau (per cent)</th><th>Inner accretionary plateau (per cent)</th></tr></thead><tbody><tr><td>0.5</td><td>43</td><td>31 (-12)</td><td>27 (-16)</td><td>25 (-18)</td></tr><tr><td>1.0</td><td>17</td><td>16 (-1)</td><td>13 (-4)</td><td>11 (-6)</td></tr><tr><td>1.5</td><td>10</td><td>9 (-1)</td><td>6 (-4)</td><td>5 (-5)</td></tr></tbody></table>

at a rate of 35–45 m Ma⁻¹ over oceanic plate offshore NW Sumatra, accreting against a continental backstop forming a wide, spectacular accretionary prism with an average fault spacing of 4–10 km. Since sedimentary layers are stratified subhorizontally in the trench, the maximum principle stress is near vertical (Fig. 9) and varies with depth, whereas the minimum principle stress is oriented almost horizontally (Tsuji et al. 2008). Therefore, the uppermost section of the trench deposits remains uncompacted and un lithified leading to a low-velocity and high-porosity zone. On the other hand, the lithification and compaction increase with depth lead- ing to higher velocities and low porosities. The orientation of the principal stress starts to become oblique towards the subduction front that increases the compaction factor horizontally, requiring the reduction in thickness of layer L1 and reverse dips of the ve- locity (porosity) contours as compared to the sedimentary layers trenchwards.

The scenario changes northeast of the trench. The direction of the maximum principle stress becomes subhorizontal to horizontal (Tsuji et al. 2008), and its magnitude increases gradually land- wards away from the trench, causing severe horizontal compaction (Fig. 10). This fact is justified by the degree of folding of the sed- imentary layers on the prism, incremental growth of ridges inside the ponded sedimentary basins (Fisher et al. 2007) and lateral de- crease in porosity on the accretionary wedge increasing the yield strength of the accumulated materials landwards. Faults associated with the accretionary prism are mostly bivergent and active in nature (since faults reach to the seafloor) and become seaward vergent (or landward-dipping) northeastwards sustaining their activity. Biver- gent faulting at the prism are also found south of our study area near the epicentre of the 2004 December 26 earthquake area (Singh et al. 2008), near the Pagai Island (Singh et al. 2011b) and in some other subduction systems, for example, Oregon and Aleutian mar- gins (MacKay et al. 1992). Development of such conjugate faulting at the subduction front is dependent on several factors. Normally high shear stress at the basal decollement and the presence of weak trench sedimentary deposits lead to form a wide taper angle, severe folding in the trench deposits and seaward-vergence (landward- dipping) thrust faulting (Byrne et al. 1988; MacKay et al. 1992). On the other hand, landward-vergence (seaward-dipping) faulting is uncommon in the subduction setting (Graindorge et al. 2008), although analogue modelling (Gutscher et al. 2001) suggests that such faulting can be developed while a thick Newtonian fluid type rheology composed of granular material overlies a ductile basal layer. Here, we imaged landward verging faults, but they seem to be linked to the bivergent nature of faulting instead, which is very pre- velant all along the Sumatra–Andaman margin (Moeremans et al. 2014).

Adjacent to our study area, Fisher et al. (2007) estimated the surface slope and the dip of basal decollement as 8°–10° and 5°, respectively, inferring the presence of a wide taper angle. Our high-resolution depth migrated images (Figs 8 and 10) show relatively

low surface slope of 6° and dip angle of the basal decollement of 4° below the frontal slope. This confirms the prism taper as very high (~10°), and consequently explains the presence of landward-dipping frontal thrust and other thrusts below the R3 and R4 ridges, wherein most of them seem to be active. Presence of such active landward-dipping faults with the high prism taper explains the high basal friction beneath the frontal slope. In addition, the velocity–depth model (Fig. 4) also illustrates the sedimentary deposits over- lying the subducting plate at the toe of the accretionary prism as high velocity (4.8–5.5 km s⁻¹), indicating their high yield strength and strong coupling with the downgoing plate, which in turn develops a suitable regime for rupture to propagate from the seismogenic zone to the subduction front and further seawards. The high-resolution image (Fig. 8) further shows that the frontal folded and faulted ridges (R1–R4) with steep slopes are clearly distinct from a wide (nearly 130 km), relatively flat accretionary plateau that lies 1600-m water depth and are separated by negative slope in the bathymetry (between 100 and 110 km offset ranges), requiring uplift of the frontal ridges or subsidence of the plateau. At the base of this neg- ative slope, there are ~3-km thick undisturbed sediments on the accretionary plateau, which seem to belong to the turbiditic sedi- ments observed in the trench, suggesting that these sediments were placed here by some process. The low taper angle at the accretionary plateau (5°–7°) implies that the basal friction beneath the S1 and S2 slopes may be comparatively less. Singh et al. (2012) correlated this drastic change in dip from slope deposits to plateau with the down- going slab geometry, suggesting that underplating of the oceanic crust and mantle megathrust as cause of differential uplift of the ac- cretionary plateau. Although the depth image (Fig. 8) does not reveal the exact plate geometry at depth, we speculate that underplating or duplexing of the oceanic plate may generate such irregularities at the plate interface. Other possibility is a 60-km wide fracture zone (F6) might be subducting in this area (Singh et al. 2011a; Carton et al. 2014) with thicker crust as compared to that fur- ther east, causing localized uplift. However, elongated high gravity anomaly all along the subduction front cannot be explained by the presence of the fracture zone or thick crust. To account for the spec- tacular width of the accretionary plateau, Fisher et al. (2007) and Gulick et al. (2011) suggested the presence of a strong material be- neath the accretionary wedge. Our high-resolution velocity–depth model (Figs 4 and 8c) and the depth image (Fig. 8) are in a good agreement with this argument. We observe surprisingly high ve- locities (4.2–5.5 km s⁻¹) at a depth of 3 km below the accretionary plateau and the depth image is chaotic due to lack of well-preserved sedimentary layers. Accreted sediments associated with the Hi- malayan orogenesis and prolonged transportation cause trench infill, bearing higher concentrations of sand and silt, (Gulick et al. 2011) which, furthermore, become more compacted due to the increase in lateral compression and dewatering, and eventually strengthen- ing the accretionary plateau. However, the size of the accretionary plateau is not constant, even it has not developed in many parts
of Sumatra subduction system. In the south, this plateau is absent (Kopp 2001; Singh et al. 2011b) due to thin sediment, whereas in north of the 2004 rupture area such marginal plateau is not well developed due to absence of suitable backstop (Gulick et al. 2011). In other areas, this plateau has been replaced by forearc high, except around 10°N between Andaman and Nicobar Islands where the forearc high has subsided (Singh et al. 2013), leading to negative slope, which could be due to the subduction on Ninety-East Ridge.

The high-amplitude reflector observed between 7.5 and 10 km depth ranges inside the trench deposits is a debated feature. On the one hand, based on previous results from south of WG2 profile, this reflector is suggested as a protodecollement and proposed to be associated with the segmentation of earthquake ruptures and tsunami generation (Dean et al. 2010), while on the other, Singh et al. (2011a) suggested that the high-amplitude reflector is a lithological boundary between pelagic sediment and turbidite. In general, development of decollements/proto-decollements are described by porosity contrast between the accreted and underthrust materials (Bangs et al. 2010), presence of weak and fluid-rich materials in the dilatant fault zones (Davis et al. 1983; Shipley et al. 1994; Bangs et al. 1999); and these horizons are identified in seismic sections by comparing the polarity reversal attributes with the seafloor reflection (Westbrook et al. 1988). To generate significant porosity contrast between the accreted and underthrust deposits, the fluid pressure should be high enough at the basal decollement. From the depth section (Fig. 9), we found that the high-amplitude layer in the trench deposits thickens landwards and it does not show strong negative polarities instead it shows variable polarities; also it is not a continuous reflector rather it is disrupted and irregular with a thickness of roughly 300 m; and most probably affected by the geometry of the top of the oceanic plate. The velocity model show high velocity at the base of the high-amplitude reflector, which gradually increases until the top of the basement without showing any inversion in velocity across the high-amplitude reflector (Fig. 11), but this could be due to the poor resolution at this depth, although Dean et al. (2010) find a low-velocity layer below this reflector using only 2.4-km-long streamer instead of a 12-km-long streamer used here. Moreover, the depth section (Fig. 8) displays that the taper is high beneath the frontal slope, which in turn, implies that basal fluid pressure may be low, resisting to develop a strong decollement surface in trench. Comparing results from further south of this study area near Nias Island, the lowermost unit overlying the top of the subducting plate was described as pelagic and turbidites (Moore & Curray 1980). Analysing other sets of data, it is observed that the high-amplitude reflectors can be found in the base of the channel-levees systems or inside fan deposits. Such high-amplitude reflectors are observed as stratified deposits in many fan deposits such as in the Bengal fan (Schwenk et al. 2005) and Amazon fan (Lopez 2001), inferring its ubiquitousness. Similar features are observed in a strike parallel profile as well (Singh et al. 2011a) and well further south (Singh et al. 2011b). Singh et al. (2011a) has observed similar very strong negative polarity reflector above that oceanic plate that has been faulted by thrust within the oceanic crust. Although we cannot rule out that this lithologic boundary could be host a decollement surface and transfer stress to the oceanic plate, the presence of low-porosity front 40 km over the oceanic plate from the subduction front suggests the locking extends up to the trench, and possibly beyond the trench over the oceanic plate as suggested by Singh et al. 2008, 2011b. The presence of steeply dipping thrust events near the subduction front further confirms this hypothesis.

# 6 CONCLUSIONS

Compiling the bathymetric observations and multichannel seismic reflection and refraction results we make the following conclusions:

(1) An unconsolidated sedimentary deposit (L1) overlies a semi-consolidated to consolidated layer (L2) that follows the underlying highly compacted to crystalline deposits (L3). The thickness of the unconsolidated layer is maximum in the trench and it decreases landwards. Velocity at the toe of the wedge ranges from 4.8 to 5.5 km s⁻¹ suggesting the sediments are highly compacted and lithified and therefore, produces strong locking favouring the rupture propagation to the front and further seawards.

(2) The HARP is irregular, disrupted and shows substantial contrast in impedance suggesting that it is a lithologic boundary.

(3) Horizontal stress is dominant in the accretionary plateau and influences the porosity distribution on the accretionary plateau. Calculations based on DEM theory for disordered cracks estimate the reduction in average porosity at the accretionary wedge with respect to the trench deposits. The porosities decrease at 0.5-km depth on average between trench and front by 12 per cent and thereafter it decreases by more 6 per cent in the inner accretionary plateau. The porosity decrease at depth in the trench suggests the vertical compaction due to high sediment thickness, whereas landward decrease of porosity is related to the high lateral compaction.

# ACKNOWLEDGEMENTS

We are highly grateful to Western Geco for funding the acquisition of a seismic reflection survey and processing of the data. We are grateful to A. Harding, N. Bangs, H. Carton, A. Chauhan, A. Arnulf and N. Hananto for their technical support.

# REFERENCES

Arnulf, A., Harding, A.J., Singh, S.C., Kent, G.M. & Crawford, W., 2011. Strong seismic heterogeneity in layer 2A near hydrothermal vents at the mid-Atlantic ridge, *Geophys. Res. Lett.*, **38**, L13320, doi:10.1029/2011GL047753.

Bangs, N.L., Shipley, T.H. & Moore, G.F., 1996. Elevated fluid pressure and fault zone dilation inferred from seismic models of the northern Barbados Ridge decollement, *J. geophys. Res.*, **101**, 627–642.

Bangs, N.L., Shipley, T.H., Moore, J.C. & Moore, G.F., 1999. Fluid accumulation and channeling along the northern Barbados Ridge decollement thrust, *J. geophys. Res.*, **104**, 20 399–20 414.

Bangs, N.L., Hornbach, M.J., Moore, G.F. & Park, J.-O., 2010. Massive methane release triggered by seafloor erosion offshore southwestern Japan, *Geology*, **38**, 1019–1022.

Berge, P.A., Fryer, G.J. & Roy, H.W., 1992. Velocity-porosity relationship in the upper Oceanic crust: theoretical considerations, *J. geophys. Res.*, **97**, 15 239–15 254.

Berryhill, J.R., 1979. Wave equation datuming, *Geophysics*, **44**, 1329–1344.

Berryman, J.G., 1980. Long-wavelength propagation in composite elastic media. 1: spherical inclusions, *J. acoust. Soc. Am.*, **68**, 1809–1819.

Berryman, J.G., 2002. Estimating rock porosity and fluid saturation using only seismic velocities, *Geophysics*, **67**(2), 391–404.

Byrne, D.E., Davis, D.M. & Sykes, L.R., 1988. Loci and maximum size of thrust earthquakes and the mechanics of the shallow region of subduction zones, *Tectonics*, **7**(4), 833–857.

Calvert, A.J., 2004. Seismic reflection imaging of two megathrust shear zones in the northern Cascadia subduction zone, *Nature*, **428**, 164–167.

Calvert, A.J., Ramachandran, K., Kao, H. & Fisher, M.A., 2006. Local thickening of the Cascadia forearc crust and the origin of seismic reflectors in the uppermost mantle, *Tectonophysics*, **420**, 175–188.
Carton, H., Singh, S.C., Hananto, N.D., Martin, J., Djajadihardja, Y.S., Udrekh, Franke, D. & Gaedicke, C., 2014. Deep seismic reflection images of the Wharton Basin oceanic crust and uppermost mantle offshore Northern Sumatra: relation with active and past deformation, *J. geophys. Res.*, **119**, doi:10.1002/2013JB010291.

Chauhan, A.P.S., 2010. Structure of the Northern Sumatra subduction megathrust using seismic reflection and refraction data, *PhD thesis*, Institut de Physique de Globe du Paris.

Chauhan, A.P.S. et al., 2009. Seismic imaging of forearc backthrusts at northern Sumatra subduction zone, *Geophys. J. Int.*, **179**, 1772–1780.

Christensen, N., 1984. Pore pressure and oceanic crustal seismic structure, *Geophys. J. R. astr. Soc.*, **79**, 411–423.

Claerbout, J., 1976. *Fundamentals of Geophysical Data Processing*. Available at: http://sepwww.stanford.edu/sep/prof/.

Curray, J.R. & Moore, D.G., 1971. Growth of the Bengal Deep-Sea Fan and denudation in the Himalayas, *Bull. geol. Soc. Am.*, **82**, 563–572.

Curray, J.R., Emmel, F.J. & Moore, D.G., 2003. The Bengal fan: morphology, geometry, stratigraphy, history and processes, *Mar. Petrol. Geol.*, **19**, 1191–1223.

Dahlen, F.A., 1990. Critical taper model of fold-and-thrust belts and accretionary wedges, *Annu. Rev. Earth planet. Soc.*, **18**, 55–99.

Dahlen, F.A., Suppe, J. & Davis, D.M., 1984. Mechanics of fold-and-thrust belts and accretionary wedges, cohesive coulomb theory, *J. geophys. Res.*, **89**, 10 087–10 101.

Davis, D., Suppe, J. & Dahlen, F.A., 1983. Mechanics of fold-and-thrust belts and accretionary wedges, *J. geophys. Res.*, **88**, 1153–1172.

Dean, S.M. et al. 2010. Contrasting decollement and prism properties over the Sumatra 2004–2005 earth-quake rupture boundary, *Science*, **329**, 207–210.

Deplus, C. et al., 1998. Direct evidence of active deformation in the eastern Indian oceanic plate, *Geology*, **26**, 131–134.

Fisher, D., Mosher, D.C., Austin, J.A. Jr., Gulick, S.P.S., Masterlark, T. & Moran, K., 2007. Active deformation across the Sumatran forearc over the December 2004 Mw9.2 rupture, *Geology*, **35**, 99–102.

Franke, D. et al., 2008. The great Sumatra-Andaman earthquakes—imaging the boundary between the ruptures of the great 2004 and 2005 earthquakes, *Earth planet. Sci. Lett.*, **269**, 119–130.

Gamage, K., Screaton, E., Bekins, B. & Aiello, I., 2011. Permeability-porosity relationships of subduction zone sediments, *Mar. Geol.*, **279**, 19–36.

Gazdag, J., 1978. Wave equation migration with the phase shift method, *Geophysics*, **43**, 1342–1351.

Ghosal, D., Singh, S.C., Chauhan, A.P.S. & Hananto, N., 2012. New insights on the offshore extension of the Great Sumatran fault, NW Sumatra, from marine geophysical studies, *Geochem. Geophys. Geosyst.*, **13**, doi:10.1029/2012GC004122.

Goldfinger, C. & McNeill, L.C., 2006. Sumatra and Cascadia: parallels explored, *EOS, Trans. Am. geophys. Un.*, **87**(52), Fall Meet. Suppl., Abstract U44A-06.

Gordon, R.G., DeMets, C. & Royer, J.-Y., 1998. Evidence for long-term diffuse deformation in the equatorial Indian Ocean, *Nature*, **395**, 370–374.

Graindorge et al. 2008. Impact of lower plate structure on the upper plate deformation at the NW Sumatra convergent margin from seafloor morphology, *Earth planet. Sci. Lett.*, **275**, 201–210.

Gulick, S.P.S. et al., 2011. Updip rupture of the 2004 Sumatra earthquake extended by thick indurated sediments, *Nature Geosci.*, **4**, 453–456.

Gutscher, M.-A., Klaeschen, D., Flueh, E. & Malavieille, J., 2001. Non-Coulomb wedges, wrong way thrusting, and natural hazards in Cascadia, *Geology*, **29**, 379–382.

Harding, A.J., Kent, G., Blackman, D.K., Singh, S.C. & Cannales, J-P., 2007. A new method for MCS refraction data analysis of the uppermost section at a Mid-Atlantic Ridge core complex, *EOS, Trans. Am. Geophys. Un.*, **88**(52), Fall Meet. Suppl., Abstract S12A-03.

Hashin, Z. & Shtrikman, S., 1963. A variational approach to the theory of the elastic behaviour of multi-phase materials, *J. Mech. Phys. Solids*, **11**, 127–140.

Henstock, T.J., McNeill, L. & Tappin, D.R., 2006. Seafloor morphology of the Sumatran subduction zone: surface rupture during megathrust earthquakes? *Geology*, **34**(6), 485–488.

Hill, R., 1952. The elastic behaviour of a crystalline aggregate, *Proc. Phys. Soc.*, **A65**, 349–354.

Hill, E.M. et al., 2012. The 2010 Mw 7.8 Mentawai earthquake: very shallow source of a rare tsunami earthquake determined from tsunami field survey and near-field GPS data, *J. geophys. Res.*, **117**, B06402, 10.1029/2012JB009159.

Huyghe, P., Foata, M., Deville, E. & Mascle, G., 2004. Channel profiles through the active thrust front of the southern Barbados prism, *Geology*, **32**, 429–432.

Jakobsen, M., Hudson, J., Minshull, T. & Singh, S.C., 2000. Elastic properties of hydrate-bearing sediments using effective medium theory, *J. geophys. Res.*, **105**(B1), 561–577.

Kopp, H., 2001. Crustal structure along central Sunda Margin, Indonesia, *PhD thesis*, University of Kiel.

Liu, C.S., Murray, J.R. & McDonald, J.M., 1983. New constraints on the tectonic evolution of the eastern Indian Ocean, *Earth planet. Sci. Lett.*, **65**, 331–342.

Lopez, M., 2001. Architecture and depositional pattern of the Quaternary deep-sea fan of the Amazon, *Mar. Petrol. Geol.*, **18**, 479–486.

MacKay, M.E., Moore, G.F., Cochrane, G.R., Moore, J.C. & Julm, L.D., 1992. Landward vergence and oblique structural trends in the Oregon margin accretionary prism: implications and effect on fluid flow, *Earth planet. Sci. Lett.*, **109**, 477–491.

Martin, J., Ozbek, A., Combee, L., Lunde, N., Bittleston, S. & Kragh, E., 2000. Acquisition of marine point receiver seismic data with a towed streamer, SEG, Expanded abstract.

Matson, R. & Moore, G.F., 1992. Structural controls on forearc basin subsidence in the central Sumatra forearc basin, in *Geology and Geophysics of Continental Margins*, AAPG Memoir 53, pp. 157–181, American Association of Petroleum Geologists.

Moeremans, R., Singh, S.C., Mukti, M., McArdle, J. & Johansen, K., 2014. Seismic images of structural variations along the deformation front of the Andaman-Sumatra Subduction zone: implications for rupture propagation and tsunamigenesis, *Earth and planet. Sci. Lett.*, **386**, 75–85.

Moore, G.F. & Murray, J.R., 1980. Structure of the Sunda Trench lower slope off Sumatra from multichannel seismic reflection data, *Mar. geophys. Res.*, **4**, 319–340.

Moser, T.J., 1991. Shortest path calculation of seismic rays, *Geophysics*, **56**, 59–67.

Mosher, D.C., Austin, J.A., Jr., Fisher, D. & Gulick, S.P.S., 2008. Deformation of the northern Sumatra accretionary prism from high-resolution seismic reflection profiles and ROV observations, **252**, 89–99.

Norris, A.N., 1985. Differential scheme for the effective moduli of composites, *Mech. Mater.*, **4**, 1–16.

O'Connell, R.J. & Budainsky, B., 1974. Seismic velocities in dry and saturated cracks, *J. geophys. Res.*, **79**(35), 5412–5426.

Saffer, D.M. & Bekins, B.A., 2002. Hydrologic control on the morphology and mechanics of the accretionary wedge, *Geology*, **30**(3), 271–274.

Scales, J.A., 1987. Tomographic inversion via the conjugate inversion method, *Geophysics*, **52**, 179–185.

Schwenk, T., Spiesz, V., Breitzke, M. & Hübscher, C., 2005. The architecture and evolution of the Middle Bengal Fan in vicinity of the active channel-levee system imaged by high-resolution seismic data, *Mar. Petrol. Geol.*, **22**(5), 637–656.

Seher, T., Singh, S.C., Crawford, W.C. & Escartin, J., 2010. Upper crustal velocity structure beneath the cen-tral Lucky Strike segment from seismic refraction measurements, *Geochem. Geophys. Geosyst.*, **11**(5), doi:10.1029/2009GC002894.

Sheriff, R.E. & Geldart, L.P., 1982. *Exploration Seismology*, Vol. 1, Cambridge Univ. Press, 253 pp.

Shipley, T.H., Moore, G.F., Bangs, N.L., Moore, J.C. & Stoffa, P.L., 1994. Seismically inferred dilatancy distribution, northern Barbados Ridge décollement: implications for fluid migration and fault strength, *Geology*, **22**, 411–414.
Singh, S.C. et al., 2008. Seismic evidence for broken oceanic crust in the 2004 Sumatra earthquake epicentral region, *Nature Geosci.*, **1**, 777-781.

Singh, S.C., Carton, H., Chauhan, A.P. S., Androvandi, S., Davaille, A., Dyment, J., Cannat, M. & Hananto, N., 2011a. Extremely thin crust in the Indian ocean possibly resulting from the Plume-Ridge interaction, *Geophys. J. Int.*, **184**, 29-42.

Singh, S.C., Hananto, N., Mukti, M., Permana, H., Djajadiharja, Y. & Harjono, H., 2011b. Seismic images of the megathrust rupture during the 25th October 2010 Pagai earthquake, SW Sumatra: frontal rupture and large tsunami, *Geophys. Res. Lett.*, **38**, L16313, doi:10.1029/2011GL048935.

Singh, S.C., Chauhan, A.P.S., Calvert, A.J., Hananto, N.D., Ghosal, D., Rai, A. & Carton, H., 2012. Seismic evidence of bending and unbending of subducting oceanic crust and the presence of mantle megathrust in the 2004 Great Sumatra earthquake rupture zone, *Earth planet. Sci. Lett.*, **321-322**, 166-176.

Singh, S.C., Moeremans, R., McArdle, J. & Johansen, K., 2013. Seismic image of the sliver strike-slip fault and back thrust in the Andaman-Nicobar region, *J. geophys. Res.*, **118**, doi:10.1002/jgrb.50378.

Taylor, M. & Singh, S., 2002. Composition and microstructure of magma bodies from effective medium theory, *Geophys. J. Int.*, **149**, 15-21.

Tsuji, T., Tokuyama, H., Pisani, P.C. & Moore, G., 2008. Effective stress and pore pressure in the Nankai accretionary prism off the Muroto Peninsula, southwestern Japan, *J. geophys. Res.*, **113**, B11401, doi:10.1029/2007JB005002.

Van Avendonk, H.J.A., Harding, A.J., Orcutt, J.A. & McClain, J.S., 1998. A two dimensional tomographic study of the Clipperton transform fault, *J. geophys. Res.*, **103**, 17 885-17 899.

Van Avendonk, H.J.A., Harding, A.J., Orcutt, J.A. & McClain, J.S., 2001. Contrast in crustal structure in across the Clipperton transform fault from travel time tomography, *J. geophys. Res.*, **106**, 10 961-10 981.

Van Avendonk, H.J.A., Shillington, D., Holbrook, W.S. & Hornbach, M., 2004. Infering crustal structure in the Aleutian island arc from a sparse wide-angle seismic data set, *Geochem. Geophys. Geosyst.*, **5**(8), 1527-2087.

Von Huene, R. & Scholl, D.W., 1991. Observations at convergent margins concerning sediment subduction, subduction erosion and growth of continental crust, *Rev. Geophys.*, **29**, 279-316.

Wang, K. & Hu, Y., 2006. Accretionary prisms in subduction earthquake cycles: the theory of dynamic Coulomb wedge, *J. geophys. Res.*, **111**, doi:10.1029/2005JB004094.

Westbrook, G.K., Ladd, J.W., Buhl, P., Bangs, N. & Tiley, G.J., 1988. Cross section of an accretionary wedge: Barbados Ridge complex, *Geology*, **16**, 631-635.

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article:

**Figure S1.** Schematic diagram illustrates the governing principle for downward continuation in time domain. The extrapolated wavefield can be calculated at any point on an arbitrary datum convolving the recorded surface wavefields with the time lagged delta functions and summing the convolved results over all the receivers for all shot gathers.

**Figure S2.** Resolution analysis is shown through the checker board test. We have created three different kinds of anomaly patterns which are different from each other in sizes both in spatial and temporal direction. The wavelengths between two anomalies decrease from pattern 1 to 3. The pattern 1 consists of the largest wavelength of 1500 m in both vertical and horizontal directions, which is a factor of 30 of our node spacing. Consequently, the pattern 3 contains the smallest wavelength of 500 m (in both horizontal and vertical directions), which is 10 times larger than our grid spacing and pattern 2 is the intermediate one that is consisted of anomalies of 1000 m. We carry out the synthetic tests in three major steps: (1) we perturb our best-fitting model (Fig. 5b) with this different patterns and the forward modelling (ray tracing) is carried out estimating the synthetic traveltime curves, which are used as the input (equivalent to observed traveltime curves) for the inversion; (2) We consider a common starting model (e.g. the original best-fitting model) for all of these patterns and the inversion is conducted and (3) When the chi-square value reduces down to 1, we stop the inversion and consider it as the final model. We subtract it from the background model (or starting model) and recover the anomalies. These three steps are followed for recovering all size of anomaly patterns. We have carried out the synthetic tests for the partly on the trench deposits and the outer accretionary wedge. The left panels of the Fig. S2 show the perturbed velocity models, which are treated as the true models and the recovered models after the inversion are displayed in the right panels. Comparing the true model with the recovered model one can estimate the percentage of reconstructed anomalies and sensitivity of the simulation. Models, which are perturbed with 1500 m anomalies, are well recovered both in shallow and deeper parts, whereas with decrease in size the degree of recovery reduces in depth but the anomalies at the shallow sections can be resolved for accretionary wedge for 500 m anomalies. From these tests, one can say that the amount of resolvability increases with depth with increasing the size of the anomalies.

**Figure S3.** The velocity-depth model derived from the high-resolution traveltime tomography of downward continued gathers is combined with the OBS tomography result using the linear interpolation. The shallow part is contributed from the high-resolution tomography while the deeper parts are obtained from OBS tomography results. This combined model is used to image the downgoing subducting plate.

**Figure S4.** (a) Input Common depth point (CDP) gather for PreSDM; (b) Output Common image point (CIG) gather after PreSDM.
(http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggu182/-/DC1)

Please note: Oxford University Press are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.