Reflection tomography by depth warping: A case study across the Java trench

. Accurate subsurface velocity models are crucial for geological interpretations based on seismic depth images. Seismic reflection tomography is an effective iterative method to update and refine a preliminary velocity model for depth 10 imaging. Based on residual move-out analysis of reflectors in common image point gathers an update of the velocity is estimated by a ray-based tomography. To stabilize the tomography, several preconditioning strategies exist. Most critical is the estimation of the depth error to account for the residual move-out of the reflector in the common image point gathers. Because the depth errors for many closely spaced image gathers must be picked, manual picking is extremely time-consuming, human biased, and not reproducible. Data-driven picking algorithms based on coherence or semblance analysis are widely 15 used for hyperbolic or linear events. However, for complex-shaped depth events, pure data-driven picking is difficult. To overcome this, the warping method named Non-Rigid Matching is used to estimate a depth error displacement field. Warping is used, e.g., to merge photographic images or to match two seismic images from time-lapse data. By matching a common image point gather against its duplicate that has been shifted by one offset position, a locally smooth-shaped displacement field is calculated for each data sample by gather matching. Depending on the complexity of the subsurface, sample tracking 20 through the displacement field along predefined horizons or on a simple regular grid yields discrete depth error values for the tomography. The application to a multi-channel seismic line across the Sunda subduction zone offshore Lombok island, Indonesia, illustrates the approach and documents the advantages of the method to estimate a detailed velocity structure in a complex tectonic regime. By incorporating the warping scheme into the reflection tomography, we demonstrate an increase in the velocity resolution and precision by improving the data-driven accuracy of depth error picks with arbitrary shapes. This


Introduction
Reflection tomography and pre-stack depth migration of multi-channel seismic reflection (MCS) data have evolved into standard seismic data processing routines in recent decades, owing to the rapid development of CPU performance and the effective adaption of seismic data processing software.Pre-stack depth migration (PSDM) is the algorithm of choice in reflection seismology to properly image steeply dipping reflectors while accounting for non-hyperbolic move-out caused by lateral velocity variations (Yilmaz, 2001;Jones et al., 2008) and thus is applied in tectonically and structurally complex geological settings in 2-D and 3-D migration strategies (Collot et al., 2011;Han et al., 2017;Li et al., 2018;Shiraishi et al., 2019).However, the quality of subsurface imaging depends on the seismic velocity model that is used for the migration.An exact determination of the velocity field is thus crucial to retrieve an optimal subsurface image.
As is well known, PSDM may also be used as a velocity estimation tool to retrieve interval velocities by performing velocity analysis on selected locations using depth focusing error analysis (Audebert and Diet, 1990;MacKay and Abma, 1992) or hyperbolic residual move-out (RMO) correction of common image point (CIP) gathers followed by a simple vertical Dix inversion (Dix, 1955) at each location independently, to manually build a structural velocity model (Audebert et al., 1997).
Manual picking, however, is not only time-consuming for seismic processors but may also lead to a subjective interpretation bias.In contrast, the velocity model design based on reflection tomography inverts all CIP locations simultaneously to update the velocity structure and yields a more complete solution (Stork, 1992;Kosloff et al., 1996).The general procedure for reflection tomography is to go into the pre-stack migrated CIP offset-domain and to measure the hyperbolic residual move-out of the depth misalignment (also called depth error) by manual picking or by automatic slowness scanning techniques (Hardy, 2003;Claerbout, 1992).Subsequently, the depth error is inverted to velocity changes to flatten the reflector signals over the entire offset range (Jones et al., 2008;Fruehn et al., 2008).For regions with a highly variable velocity, however, a more severe non-hyperbolic depth misalignment becomes a common situation, especially with increasing source-to-detector distances.To guide the velocity changes along the subsurface structures, predefined horizons can be included during the picking procedure as preconditioning during a layer-based tomography (Riedel et al., 2019).An overview of the model building techniques can be found in Jones (2003).
Several workflows have been established for velocity estimations depending on the different acquisition types and seismic wavefields available.Long streamer acquisitions with respect to the target depth in shallow water depth offer the possibility to invert for complementary wavefields, the near-vertical reflected events and the horizontal propagating refracted arrivals from the same dataset.Gras et al. (2019) used selected reflected arrival times and the refracted arrivals for a travel time tomography to estimate an initial velocity.This velocity was used for a subsequent full wavefield inversion followed by PSDM.
To image crustal structures in a deep-water environment, Górszczyk et al. (2019) used complementary wavefields from streamer data and ocean bottom stations (OBS) recordings.A first arrival tomography of the OBS data produced an initial velocity for a full waveform inversion of the OBS recordings.To minimize the RMO of CIP gathers of the PSDM data with Deleted: on the waveform inversion velocity, a slope tomography of local reflector elements was used to further improve the velocity model for a final migrated image.
An established processing flow for reflection tomography of MCS data includes the determination of an initial P-wave velocity field, multiple attenuation with adaptive subtraction (Verschuur et al., 1992;Guitton and Verschuur, 2004), Kirchhoff PSDM, depth error estimation of CIP offset-gathers, followed by ray-based tomographic inversion to update the velocity field.
In general terms, the tomographic inversion identifies an optimal model that explains the observed input data (Bishop et al., 1985).An equivalent approach is the migration velocity analysis technique which is based on forward and reverse ray-based propagation of travel times to find a velocity model that minimizes the CIP depth error (Van Trier, 1990).While the Dix inversion strips off the layers from top to bottom in a flat-layer approach, the ray-based tomography accounts for dipping layers and lateral velocity changes within the streamer length (Jones, 2010).Both the pre-stack depth migration and depth tomography algorithms rely on the definition of the initial seismic velocity field in the subsurface (i.e. the starting model).
By relating changes of the CIP depth errors to changes in velocities along source-receiver rays in the direction normal to the local reflector dip, a new velocity can be calculated to minimize the CIP depth errors.To solve the general non-linear inverse tomographic problem, the velocity error is gauged iteratively, inverted, and updated as depicted by the loop in Fig. 1.

Multi-Channel Seismic Data
To circumnavigate these issues and increase the picking accuracy, we applied a warping technique called "Non-Rigid Matching" (NRM).By calculating the depth error shift of seismic trace samples by comparing neighbouring traces along the complete offset of closely spaced CIP gathers, we improve the depth error estimation without any hyperbolic assumption or predefined depth horizons of the subsurface structure.
Here we present the NRM technique for the depth error estimation as a pure data-driven automatic picking method.We demonstrate the advantages and limitations of the NRM method using a synthetic gather.We then apply a combination of NRM with ray-based reflection tomography to field data of pre-stack depth migrated seismic sections from the Sunda convergent margin offshore Lombok island, Indonesia.As initial velocity, a wide-angle tomographic inversion of a collocated 2-D OBS seismic line was used.The reflection profile is characterized by an accretionary prism of strongly folded sediment with limited reflector continuity, which makes manual velocity estimation extremely challenging.
2 Method: Non-rigid matching and reflection tomography

Non-rigid and warping matching techniques
The non-rigid matching (NRM) or "Warping" methods are computer-based image matching technologies that aim to estimate a flow pattern (vector displacement in 3-dimensions) of a sequence of images with additional smoothness constraints (Horn and Schunck, 1981;Wolberg, 1990).Compared to a rigid matching like translation, rotation, or even affine transformation, NRM is developed to handle situations when the transformation is non-linear (Pappu et al., 1996).The benefit of NRM regarding the non-linear transformation substantially improves seismic imaging and inversion methods through matching and tracking horizontal and vertical displacements of seismic events with high accuracy in the depth and time domains.
NRM or warping applications were first introduced for 3D time-lapse seismic data by comparing two seismic cubes acquired at different acquisition times with a special focus on depth formation changes resulting from hydrocarbon production (e.g., Rickett and Lumley, 2001;Aarre, 2008;Tomar et al., 2017).The image displacement warping method of Hall (2006) estimates a full 3D local vector deviation employing an iterative search of maximum correlation using a deformable mesh for sensitivity and quality analysis, whereas Hale (2009Hale ( , 2013) ) based his dynamic image warping (DIW) on 1D cross-correlation optimization schemes in each dimension to estimate the vector displacements.By solving a set of 1D equations and separately including spectral whitening and a Gaussian low-pass filter, a stable 3D solution is achieved iteratively by minimizing the difference of the reference and the current wavefield corrected by the estimated displacements.This method is able to calculate rapid and large shifts, both in time/depth and in space, and overcomes the restrictions of limited shifts due to time/depth windowing used by cross-correlation methods (Zhang et al., 2014).In contrast, the NRM method introduced by Nickel and Sønneland (1999) uses 1D Taylor expansions for each vector component, which are separately solved for each dimension to converge to a 3D solution by minimizing the difference of the reference section and the current wavefield corrected by the estimated displacements like by the warping method of Hale (2009) (Aarre, 2006;Aarre, 2008).To stabilize the results, additional Deleted: Good quality of imaging and of the tomographic result must fulfil both the conditions of accuracy and precision of the input 125 information, which describe the trueness and density of the input depth picks of the depth error (Jones, 2003).The precision may be improved by setting a smaller vertical and lateral picking interval to maximize the reliability of the tomography.On the other hand, accuracy is strongly limited by signal interference, background noise, 130 side echoes, and accurate depth error information.To circumnavigate these issues and increase the picking accuracy and precision, we applied an innovative constraints are implemented, e.g., bandlimited application, smoothing for spatial continuity, and avoiding vertical shifts that would swap neighbouring depth samples.

Deleted
A number of new geophysical applications for pre-stack event tracking by the warping technique have been introduced in the scientific community in the last decade (e.g., Perez and Marfurt, 2008;Reiche and Berkels, 2018;Sripanich et al., 2020).
The main objective of all these applications is to efficiently define a reference data ensemble and calculate the displacement 145 shift from any data ensemble to match the reference data ensemble.The unique selection of a reference ensemble depends on the individual purpose of the application.Perez and Marfurt (2008) estimated vertical and spatial displacements with a modified cross-correlation method from Rickett and Lumley (2001).A displacement estimation between 3D common angle binned migrated sections to a reference stacked volume was used to improve the stack quality and resolution.Reiche and Berkels (2018) sorted the migration data into common offset sections and selected the smallest offset section as a reference 150 section, and calculated the displacement from all other offset sections to the smallest offset section in order to calculate the move-out curvature and flatten the common-mid-point (CMP) gather.Sripanich et al. ( 2020) estimated reflection move-out dip slopes on 3D CMP gather directly by a plane-wave destruction filter (Fomel, 2002) to flatten events by nonstationary filters.This processing sequence can be seen as a warping process by an application of time-variant static corrections.
The wide range of possible applications for lateral and vertical displacement estimations even in three dimensions make the 155 NRM and DIW attractive e.g.quantitative estimations of vertical and horizontal reflector shifts due to repositioning of remigrated data.As both methods are very similar by iteratively minimizing the difference of two sections, we use only the NRM for the application of estimating the RMO.The calculation of vertical shifts between a reference section to an actual section is the most simple's application and allows us to compare the results to a plane-wave destruction filter (PWD), which is commonly used for estimations of reflector dips.160 The application of NRM for the vertical displacement calculation in a CIP gather requires a reference CIP gather, similar to a time-lapse application.This can be achieved with a relative reference scheme by duplicating the current gather and shifting the traces laterally to larger offsets by one trace position to form the reference gather.The calculated relative displacement shift between the two gathers correspond then to vertical and spatial smooth event slope dips of a trace to its previous trace.
As neighbouring traces of a CIP gather show a strong similarity of the waveform and amplitude without spatial aliasing, the 165 NRM gather matching can then be used to estimate the vertical displacement also known as depth error without any physical assumptions.In this way, the application of NRM for CIP gathers overcomes the limitation of residual move-out estimations inherent in conventional semblance scanning of predefined functions like linear, parabolic, or even higher-order curvatures.

NRM synthetic data example
For field MCS data, due to the complex subsurface structure and seismic acquisition geometry, as well as the anisotropic 170 physical world, three main unique classified situations represent the main difficulties for analysing the residual move-out.To test the advantages and shortcomings of the warping method we created a synthetic CIP gather in Fig. 2a.The gather consists of three sets of events including 0.1% background noise of the maximum amplitude from top to bottom: (1) a symmetrical  et al., 2017).A recursive application of the depth variant static displacements for each trace will flatten all events of the current CIP gather Deleted: Generally, we introduce here a relative referenced scheme 185 that always compares and calculates the relative displacement shift (also known as the event slope dip) between neighbouring traces, rather than comparing all traces to a near-offset reference trace as is commonly done in the conventional approach.With the existence of noise and independent events, the relative referenced NRM rejects 190 possible outliers due to predefined smoothness constraints and determines an optimal spatial transformation for mapping one trace to the other while preserving a similar salient structure along the traces.The stability and compatibility of this application to high noise-level data are some of the main advantages of non-rigid matching 195 compared to conventional semblance scanning.¶ A number of geophysical applications for pre-stack event tracking by the warping technique have been introduced in the scientific community in the last decade (e.g., Reiche and Berkels, 2018;Sripanich et al., 2020).

200
Deleted: The main objective of all these applications is to efficiently define a reference trace or dataset and calculate the displacement shift from any trace to match the reference trace.The unique selection of the reference trace section depends on the individual purpose of the application.Reiche and Berkels (2018) 205 sorted the migration data to common offset sections and selected the nearest offset gather as a reference section, and calculated the displacement from all other offset sections to the nearest one in order to calculate the move-out curvature and flatten the common-midpoint (CMP) gather.

220
synthetic gather based on these three extreme situations (Fig. 2).
lateral shifted diffraction-like event, which is unrealistic but was included because it cannot be approximated by any linear, parabolic or hyperbolic curvature; (2) two intersecting parabolic curvature events with opposite polarity.To get the undercorrected positive polarity event horizontal aligned the velocity above of must be reduced whereas the velocity above the overcorrected event with the negative polarity must be increased.This situation occurs if the background velocity model is not 225 well adapted to the data e.g. a vertical velocity increase between the two events; and (3) one parabolic event with a local curvature anomaly, offset-dependent wavelet amplitude, and frequency attenuation.To get the under-corrected general trend of the event horizontal aligned the velocity above must be decreased.To further align the local curvature anomaly the velocity above must locally be increased to fully flatten the signals over the complete offset range.
Because the NRM displacement field in Fig. 2b is calculated in a relative referenced scheme of a trace to its previous trace, 230 the dip field contains relative dip displacements.The red colour of positive values in each trace shown in Fig. 2b suggests that a corresponding trace sample in Fig. 2a should be shifted downward to match and align to its previous trace sample.Bluecoloured negative values require shifting in the opposite (i.e., upward) direction.A zero-displacement value appearing at the apex of the symmetrical diffraction events illustrates the fact that the dipping angle at this location of the event is zero.The NRM field of these three sets of synthetic events follows the general local dip trend well.235

Depth variant alignment from relative displacement correction
Since the NRM field contains the full information of the relative depth variant shifts of the seismic events, the NRM field can be used to flatten the input seismic section, which has several advantages for the depth error calculation and as quality control of the validity of the displacement field.Intuitively the second trace of Fig. 2a ), (1) 245 where  [",:] is the original synthetic seismic trace array at the i th trace, ℎ [",:] represents the NRM displacement field at the i th trace, the index i represents the actual trace number index of the dataset ( ∈ {1, 2, 3, 4, … … }), and the function  &"' represents an irregular linear interpolation function.The function  &"' for any corrected sample  [ ",-] could be expressed as in Eq. (2): Deleted: In Fig. 2a, three sets of synthetic seismic sections are generated from top to bottom, which are first a symmetrical diffraction-like event, second two intersecting events with opposite polarity, and third a seismic event with an offset and amplitude varying wavelet and some non-linear local undulation.¶ The NRM depth variant shifts of the second trace are applied to the second trace and all following traces in the gather.This sequence repeats until the last trace is corrected.
where the m and n are the closest irregular (non-integer) index to j in the depth corrected index array  [ : ] by the NRM field for any given seismic trace  [ ",: ] .Please note that the displacement shifts from seismic amplitudes have much higher precision than the traces' depth sample rate, and the correction of the depth indexes will end up with non-integer numbers in the intermediate index array  [ : ] .Therefore, an irregular linear interpolation is needed in the calculation to recover the original regular depth index j.The depth index j will be immediately obtained after the linear interpolation.does not have to be an integer.Thus, the array  [ : ] could be simply expressed as: where the ℎ [ ",: ] represent the NRM displacement field at the i th trace.By applying the Eq.(1-3), one could readily get the flattened synthetic seismic section  [ :,: ] (Fig. 2c).The Python programming with the flattening of the seismic section in Fig. 2c is documented in the online repository.
The NRM flattened events by the recursive depth variant correction in Fig. 2c provide quality control of displacement calculation for these three special situations.Generally, the squeeze-and-stretch effect of the non-linear displacement correction is inevitable for a multi-trace gather, but as shown here, the displacement shift correction adequately dealt with most of the simulated examples.The first symmetrical diffraction events get optimally flattened, with no significant change of the wavelet shape.The third non-linear undulation with a wavelet variation effect is effectively flattened at the peak amplitude, but strong wavelet stretch is visible.After the NRM displacement correction, the wavelets at mid offsets (1.5 to 2.5 km) get squeezed, and at far offsets (3.5 to 4.0 km) stretched equivalent to normal moveout stretch effects.The crossing region of the two intersecting events is flattened well but suffers from a significant stretch effect, which introduces substantial artificial lowfrequency energy between the two events between offset 2.5 and 3.

RMO automatic picking by tracking through NRM displacement field
As the flattened NRM relative displacement field contains the information of relative depth shift from trace to trace along depth slices, the estimation of the depth error is achievable by predefining a start tracking depth at the first trace (nearest offset) and analyzing the corresponding depth slice of the flattened NRM field in Fig. 2d.Given a pre-defined starting pick  7 at the nearest offset, the change of the residual reflector depth  [:,82] , which is needed for the reflection tomography, can be extracted along the flattened NRM displacement field ℎ [:,:] by calculating a cumulative summation ℎ [:,:] along the depth slide at depth where the array  . (5) The array ℎ [",:] represent the i th flattened NRM trace.With the knowledge of the residual reflector depth  [",82] , one could readily get the RMO of any reflector in a seismic gather: where array  [:] represents the absolute depth of an RMO sequence over the gather for a series of continuous reflectors.By applying the Eq.(4-6), and using the simple synthetic seismic (Fig. 2a), the four series of RMO picks illustrated in Fig. ( 2e) represent the auto-tracked RMO depth of the events.The Python programming of the Eq.(4-6), the calculation of the RMO of a series of continuous reflectors (Fig. 2e), and the absolute NRM field calculation (Fig. 2f) are documented in the online repository (https://github.com/xyywy/Sup-SE2021-40)as web-based interactive script files.
As a result, all residual move-out (RMO) depth error picks follow the amplitude peaks of the seismic events except for the "X" shaped interfering reflectors.The NRM displacements are misled by the crossing point and switch to an event that should not be followed.This kind of "V" shape depth error information will undermine the reliability of the tomographic result because the depth error branch is a combination of two different events.To avoid this mismatch a quality factor is introduced and assigned to each individual pick along the depth error branch.A sliding trace summation window along a depth slice or by comparing a near-offset stack to the individual events can be used.In this example, a sliding trace summation would detect a rapid decrease of quality due to the polarity change resulting in a destructive summation.Picks with lower quality as a predefined threshold can be deleted.In this special case, the two depth error branches of the intersecting events will be split into four individual depth error branches which the reflection tomography will handle correctly as four independent reflection events.If these crossing events happen as a result of interfering noise like surface-related or interbed multiples, the unwanted events should be attenuated prior to NRM by a dip filter in the CIP gather (Fig. 1).
To verify the accuracy of the NRM method we compared the picking results of the NRM method to a plane-wave destruction filter PWD method (Fig. 3).The PWD method is splitting the data into spatially and temporal windows and assumes that the slopes are stationary within each window (Fomel, 2002).In contrast, the NRM method is iteratively minimizing the energy difference between an adjusted gather and a reference gather for all events simultaneously.For both methods, the estimated depth error picks in Fig. 3a are near the maximum amplitude peak of the events.On the far offset traces with reduced amplitudes the NRM show less accuracy (Fig. 3b).Due to the strategy by the NRM method to minimize the energy difference, the strongest amplitude events will dominate this inversion result.This strategy is useful by comparing two seismic cubes with a special focus on depth or spatial shifts to reduce distortions from low amplitude events.To avoid small-amplitude events in a CIPgather to be underdetermined in their dip estimation an additional gain balancing is recommended to be applied before an NRM application.

Effective RMO selection based on semblance analysis
Giving predefined starting depths for the RMO tracking is applicable in a simple synthetic test which contains only few continuous reflectors.However, it is unrealistic to define each reflector's starting depth in real data or complex synthetic examples, which could have several tens of effective reflectors in one migrated CIP gather (Fig. 4a).An efficient approach to determine the reflectors' start-tracking depth is by analysing the flattened CIP-gather (Fig. 4c) by a semblance-weighted gridbased scheme.
A Semblance s[j] value, which is a quantitative measure of the similarity of a number of traces (Yilmaz, 2001) in the seismic section, is described as: where n is the maximum number of traces,  [",-] represent the seismic section, the index i represents the trace number, and the index j represents the sample depth.
By conducting semblance calculation on the NRM-flattened seismic section (Fig. 4c), the flattened section could not only provide NRM field's quality control but also sheds light on selecting effective reflectors and determining the starting depth for RMOs' auto-tracking.By calculating the semblance value for each depth sample along the flattened seismic section, one can 410 reject unwanted picks, by setting up a threshold of semblance limitation (e.g.0.5) in depth.In the synthetic example with several reflectors (Fig. 4a), RMO picks are digitized in zones of good reflector continuity (semblance > 0.5, see coloured dots in the left panel of Fig. 4c) and RMO picks are rejected in non-reflection or weak zones (semblance < 0.5) (Fig. 4d).The math application of this auto-picking scheme is documented in the online repository (https://github.com/xyywy/Sup-SE2021-40). Because the plane wave destruction filter (PWD) is widely used to estimate the move-out dip slopes on seismic sections or gathers (Fomel, 2002;Sripanich et al., 2020), we applied the outlined processing sequence of displacement corrections based on a PWD filter in the online repository as an alternative method to the NRM method.

Methodology of the ray-based grid tomography with CIP depth errors
The basic concept of iterative ray-based grid tomography is to find subsurface velocity perturbations that minimize the residual moveout depth error picks of the initial migrated CIP gathers (Woodward et al., 2008).Conditions that must be fulfilled are preserved arrival times.The calculated arrival time  of ray path for a given depth error pick  through the initial velocity model must be preserved in the updated velocity model.That implicates that a change of the residual reflector depth  must be compensated by small changes in the velocity model  " where the index i corresponds to a grid node in a gridded model.
For an acoustic reflection it follows the residual migration equation of Stork (1992): Here t is the preserved arrival time, z is a change in reflector depth,  is half the opening angle between source and detector rays at the reflector,  is reflector dip, and  the velocity at the reflector depth of the actual model. " is a change in the velocity, and / " is a change in travel time corresponding to a change of velocity  at grid node i.As / " is calculated independent of the ray parameters  and , the ray path bending is assumed not to change during a velocity update.From this follows that only small velocity perturbation should be estimated for each iteration step.
The CIP tomography must find the velocity change  " that is needed to flatten migrated reflectors and eliminate the picked reflector depth error for each offset ℎ to a new depth  ) B based on Eq. ( 8): Due to the unknown residual migrated depth  ) B , the CIP tomography minimizes the difference ( ) B −  7 B ), where ℎ=0 is the nearest offset of a picked depth error branch, not necessarily zero offset, and ℎ a non-nearest offset.This yields to Eq. 10): We have a set of linear equations for each pick  ) at offset ℎ along a depth error branch relative to the nearest pick  7 for offset 0 along the depth error branch and that for many depth error branches along the whole profile.The tomography equation in matrix notation is written in Eq. ( 11) where  weights individual depth error areas,  is the matrix of the ray path term calculated by the residual migration term in brackets of Eq. ( 10);  is a scale length smoother with the predefined wavelength in lateral und vertical direction;  together is the velocity update vector, W is a damping factor allowing to adjusts the update magnitude of the model, and  is our accumulating picked depth errors between non-near and near-offset picks.The aim of the reflection tomography is now to solve the equation to find a  that will explain the residual depth errors .Additional regularisation and weighting schemes need to be applied to find a velocity change  that will minimize Eq. ( 11).Details can be further seen in Woodward et al., 2008.
The choice of the parameters of weight , smoothing scale length , and damping factor W are strongly data dependant.
Tomographic inversion works by iterative velocity updating to minimize the observed residual velocity error (Eq.11).In the CIP gathers, the depth errors are distributed over the entire offset range and yields the estimation of the interval velocity changes along their ray paths between the source and receiver (Jones, 2010).Due to the linearized approximation in Eq. ( 10) by ignoring the ray path bending during the estimation of the velocity update, several iterations with only small velocity updates (e.g.~10%) are recommended.Additionally, a large volume with a high spatial density of the residual depth errors needs to be picked in order to stabilize the linear equations by the redundancy of information (Jones, 2003).Once there is a conflict occurring between some of the equations in one grid cell, the minority picks (which could be good or bad) will be rejected by the tomography algorithm in order to get a stable and self-consistent result.Unrealistic picks have an unfavourable effect on the tomographic results when they become the majority.Besides the pick density and depth error accuracy , the smoothing scale length  is the most important parameter for the grid tomography.
If the initial velocity is not well determined, e.g.smoothed depth converted stacking velocities or depth converted pre-stack time migration velocities, a long wavelength to short wavelength velocity update is a preferred strategy by reducing the smoothing scale length  in Eq. ( 11).The first iteration is applied with a spatial smoothing length covering at least twice the CIP ray path coverage aperture to update the long-wavelength velocity structure only.In the following iterations, the scale length is successively reduced for each iteration to receive increasing velocity details as shown in Woodward et al., 2008.
If an initial background velocity is well determined e.g.depth focussing analysis from pre-stack depth migrated data, each iteration can be applied immediately with multiple scale lengths, starting from the longest to shortest scale length for each velocity update.Independent of the grid-based reflection tomographic inversion strategy it is common to stop the iterations if an iteration does not contribute any more to the flatness of the CIP residual moveout.Qualitative control will give a comparison of the reflectors horizontal alignment in the CIP gather with respect to a previous iteration or initial iteration.A more quantitative measure to stop an inversion is to define a lower limit of percentage velocity change which must be achieved (e.g.

3%).
In  crossing the subduction trench and accretionary wedge (yellow line) show highly complex structures where standard velocity analyses mostly fail by discontinuous highly dipping structures.
In contrast to the processing by Lüschen et al. (2011), where the velocity model was built iteratively from interpreted depth focussing analysis in a top down approach, we used a collocated wide-angle tomography model as initial velocity in combination with a data-driven grid-based reflection tomography.The use of wide-angle and refracted velocities may be 500 influenced by anisotropy, but gives the most confident velocity in the deeper subsurface due to the limited streamer length.

Study area and MCS data pre-processing
The multi-channel 2-D reflection seismic profile BGR06-313 that we use in three field examples was using a 3000 m long, 240 channel digital streamer with a group distance of 12.5 m at a towing depth of 6 m.A two string G-Gun array of 3080 in 3 (50.8l) volume with a nominal shot point distance of 50 m was used as a source across the southern Java trench in the south-505 eastern part of the Sunda subduction zone (Lüschen et al., 2011) as part of the SINDBAD project during RV SONNE Cruise SO190 (Fig. 5).
The seafloor depth ranges from 1.5 km near the shore on the northern part of the line to 6.5 km in the deep-sea trench.Details of our seismic processing sequence are provided in Tab. 1.In preparation for the Kirchhoff PSDM, the multiple reflections have been attenuated using a free surface multiple prediction (Verschuur et al., 1992) followed by a frequency-split 2D adaptive 510 least-square subtraction (Robinson and Treitel, 2000;Guitton and Verschuur, 2004), and a Radon transform dip filter (Hampson, 1986).

Initial velocity building from wide-angle tomography
The initial velocity model for the reflection depth tomography was merged from an OBS velocity tomographic inversion of a

Deleted: ocean bottom seismometers (OBS)
Deleted: manually 565 estimated from the near seafloor structure at coarsely sampled CMP locations by interactive semblance velocity analysis.This near subseafloor velocity adjustment was needed because the MCS reflection and OBS acquisitions were split into two cruise legs and both profiles did not completely coincide as seen by the mismatch of the seafloor depth at the lower slope and the trench axis (Fig. 6a).Due to a gap of three OBS positions in the trench axis the velocity structure was not well determined near the trench axis with lower slope sediment velocities of more than 3800 m/s between CDP 25000-26500 at a depth of 7000 m 570 (Fig. 6).At a first step, the MCS velocity analysis was based on pre-processed CMP gathers and interactive picked on semblance analysis with an increment of 500 CMP locations.Based on this smoothed stacking velocity a pre-stack time migration was subsequently applied and a second interactive semblance velocity analysis on the migrated CIP with the same increment delivered a smoothed and depth converted velocity model for the upper 2 km below the seafloor.To finalize the initial tomography model building, the adjusted velocity at shallow depth was merged with the wide-angle velocity model and 575 used for the following NRM tomography (Fig. 6b). the trench axis where side reflections and cross dipping structures due to the rough seafloor topography were observed (Fig. 610 5).A regional depth variant water velocity (Tab.2) was extracted from the Climatological Atlas of the World Ocean MB-System (Levitus, 1982) and used for the entire profile.
It should be noted, that the wide-angle velocity could not be updated at a depth greater than 4 km below seafloor during the initial velocity building and the subsequent reflection tomographic inversion since the limited streamer length (3 km) does not provide enough residual move-out sensitivity.615

Reflection tomography attribute data
For the automatic residual NRM picking (see Fig. 4c, 4d from the synthetic example) an initial depth slice increment of 50 m for the residual depth error pick tracking through the NRM displacement field was depth adjusted based on minimum threshold semblance values by scanning along offsets of the CIP gather.
One attribute presented for the tomography is the reflector dip field  (Eq.8) of the migrated section, which determine the 620 ray coverage propagation direction (see Fig. 7e as an example).A second attribute is the spatial variant weight function  (Eq.9) calculated from the spatial coherency of the depth migrated structure to weight the picks of the depth error branches (see Fig. 6f as an example).

Data examples
Each iteration loop in the tomographic processing flow (Fig. 1) included six sequential applied scale lengths smoothing's  625 (Eq.9).Starting from the longest down to the shortest application sequence, each smoothing was applied over the complete depth range (Tab.3).In total 5 iterations of velocity updates were applied, where for the following presentation only the initial and the final results are shown for the purpose of comparison.
In the data examples, we show three different structural settings with results of the velocity model, the corresponding PSDM sections, and the NRM displacement field, and the spatial coherence field together with the reflector dip field of the final 630 migrated section.To document the change in the CIP-gather domain in detail, we additionally compare selected subareas of initial and final CIP gathers, the calculated NRM displacement fields, the residual depth error picks, as well as on overlay display with the CIP gather and the depth error picks.

Sediment basin NRM tomography
The first field data example, "Example A," at the northern end of the profile (Fig. 7), is a shallow sediment basin with layered 635 interfaces and continuous reflectivity and represents an optimal site to obtain a reliable velocity model in a 2-D multi-channel seismic survey.A CIP-gather increment of 32 (200 m) was analysed along the profile with the NRM method.In total, five iterations of tomography loops (Fig. 1) were applied to this data example.An enlarged view of the initial velocity model  ranging from CDP 46700 to CDP 50800 is displayed in Fig. 7a.The resulting initial Kirchhoff pre-stack depth migration (Fig. 7c) retrieves a coherent image of the shallow sedimentary portion, while the energy in the deeper part close to the basement is 685 not very well collapsed, resulting in a series of over-migrated events.The displayed reflector dip field (Fig. 7e) and coherency field (Fig. 7f) are extracted from the final migration section.
The reflector dip is used for the ray propagation direction during the tomography, and the coherency field is used as an additional weighting of RMO depth error picks in spatial coherent subsurface areas.The two attribute fields were recalculated for each iteration of the tomography loops (Fig. 1).After five iterations of the NRM based depth tomography and Kirchhoff 690 PSDM, the reflection energy is much better collapsed and shows more focussed and continuous signals, especially in the deeper part between 5.2-5.6 km (Fig. 7d).Furthermore, the final velocity model (Fig. 7b) displays lateral velocity variations that mimic the form of the base of the sediment basin.This is well demonstrated by the 3000 m s -1 velocity contour that mimics the shape of the boundary between the highly reflective basement (below) and the less reflective but more laterally continuous reflections of the sedimentary sequence (above).695 Moving into the pre-stack CIP domain, a series of CIP gathers ranging from CDP 46700 to 50800 (same profile range as in Fig. 7) are selected and displayed in Fig. 8a with an increment of 32 (200 m).A dip filter is applied to the gathers to eliminate the extreme dipping events and migration noise.The NRM field in Fig. 8c shows the initial relative displacement values for each data sample.The information below the basement is muted by a digitized basement horizon.The distinct block of blue colour within the red rectangle in Fig. 8c, at a depth of 5.0 km to 5.8 km, illustrates a general velocity overestimation in the 700 overlying sediment.The RMO depth error picks calculated from the NRM displacement field, as a data-driven automatic picking method without any assumption of its curvature, is the main input information for the tomography (Fig. 8e).Figures 8b, 8d, and 8f show the final flattened CIP gather, NRM displacement field, and RMO depth error picks, respectively.
Compared to the initial data, the updated events in the CIP gather become optimally flattened.The depth of the basement shifts upwards by 0.2 km due to the velocity reduction of the final model.In the final NRM field (Fig. 8d), the velocity 705 overestimation error in the region of the red rectangle is substantially reduced.However, some residual move-out undulations from the initial to the final stage remain, as seen in detail in Fig. 8g and 8h from the CIP-gathers overlain with the RMO depth error curves.Ideally, the final NRM displacement field in Fig. 8d should have no NRM depth shift anymore, and all depth error picks should align horizontally.This cannot always be achieved, as the tomography finds only the solution that minimizes the depth error with respect to the smallest scale lengths (Tab.3).710

Accretionary wedge NRM tomography
In the following, the PSDM profile (location marked by the yellow line in Fig. 5) with the final velocity model overlain in Fig. 750 9 will be further analysed.
Two data examples (Fig. 5, examples B and C, marked in red) in the blue rectangles within the accretionary wedge in Fig. 9 show distinct levels of complexity.The selected upper slope area is characterized by strongly folded continuous reflector sequences, whereas the lower slope area contains only short reflector segments with varying dips.The lack of coherent reflective signals in this highly deformed accretionary prism leads to a severe difficulty to accurately evaluate the residual 755 move-out in the CIP gathers, especially if a constant spatial analysis increment (e.g.500 CDP = 3125 m) is greater than or equal to the lateral dimensions of velocity structures to be resolved (e.g. the piggy-back basin CDP 30500 -31000, 3 to 4 km depth in Fig. 9b).As a consequence of the spatially complex reflectivity pattern, the automated CIP analyses were reduced to an increment of 16 (100 m) to achieve more redundancy of depth error estimations during five iterations of the tomography.

Upper slope NRM tomography 760
Our second field data example focuses on a sequence of thick sediment tilted by compressive deformation in the region marked by the blue rectangle example B in Fig. 9b. Figure 10 provides a detailed image of the PSDM section and velocity model from the initial and final stages.The final velocity (Fig. 10b) is significantly reduced compared to the initial velocity model (Fig. 10a) in the shallow part and significantly increased compared to the initial model at depths of 5.2-6.2km.The reflector sequences of the anticline structure between CDP 29300 and 29500, from 4.0-4.4km depth, are more continuous in the final 765 image (Fig. 10d) than in the initial image (Fig. 10c), especially at the top of the anticline.The dip of the folded reflector sequence between CDP 29800 and 30100, above 4.8 km, is more continues dipping in the final image (Fig. 10d), since the residual depth error is better flattened (Fig. 11g and h), and the reflector dip in the PSDM section increases steadily with increasing distance from the apex of the fold (Fig. 10d).By contrast, the initial image in this same region (Fig. 10c) shows an abrupt change in the dip near the apex of the fold.770 Comparing the initial and final CIP gathers in Fig. 11a and 11b inside the red rectangle, strong downward dipping reflections indicate the requirement to reduce the initial velocity significantly.The NRM displacement field in Fig. 11c provides a more quantitative view of this requirement, seen by the strong blue colour with more than 2 m depth error per trace distance.The calculated RMO picks in Fig. 11e and overlain on the seismic image (Fig. 11g) follow the seismic down dipping reflection trend quite accurately.After the tomography, the final NRM displacement is significantly reduced (Fig. 11d and 11f), and the 775 residual calculated depth error in the red box (Fig. 11h) is reduced, and the reflectors are in better horizontal alignment.To the left of the red box between CDP 29600 and 29800, above 4.4 km depth, the tomography could only partially remove the depth error (compare Fig. 11c and 11d).The reflections in this region could only be aligned with velocities far below the 835 water velocity, indicating that side echoes or cross dipping structures in this region prevent a reliable subsurface velocity determination.To avoid such unrealistic velocity updates during the tomography, a minimum velocity of 1750 m s -1 below the seafloor was defined as a precondition.

Lower slope NRM tomography
In the lower slope region (Fig. 9a, example C), sediment layers are segmented and folded as a result of the regional compressive 840 deformation exerted by the subduction accretion processes.The initial pre-stack depth migration example is shown in Fig. 12c.
After the tomography, the final velocity increased by 500 m s -1 on average (Fig. 12b), resulting in a significant increase of the velocity gradient compared to the initial velocity model (Fig. 12a).In the final PSDM section (Fig. 12d), the reflector strength generally increased, and new reflector segments became emphasized compared to the initial migration (Fig. 12c).This is especially evident in the depth range from 6.0 to 6.8 km.845 In the initial CIP gathers displayed in Fig. 13a, the reflector distribution appears largely uncorrelated, and no clear trends are visible, particularly within the red box.In the initial NRM displacement field (Fig. 13c), there is a general positive depth error character that dominates the gathers, as indicated by the red colour, especially within the red rectangle and in the initial residual depth error (Fig. 13e).
By increasing the velocities based on the tomography result, this misalignment is reduced both in the final NRM displacement 850 field (Fig. 13d) and in the final residual depth error illustrated by the generally more horizontal alignment of the events (Fig. 13f).In the enlarged view of Fig. 13g and 13h, the general positive dip trend has been mostly removed.However, local reflector misalignment is still observed, as documented by the local blue colour in the NRM displacement of downward dipping events (Fig. 13d).Even after the tomography, the two local anomalies of four neighbouring CIP gather between CDP 28200 and 28400, at 5 km and 5.6 km depth, were not correctly aligned.These local anomalies have a lateral dimension of ~200 m 855 and are therefore three times shorter than the smallest horizontal scale length smoothing used for the last iteration of the tomography (Tab.3).

Final velocity model and reflectivity structure 950
The final depth image with the velocity model (Fig. 9) and the final subsurface velocity model (Fig. 14b) compared to the smoothed initially velocity (Fig. 14a) show local velocity changes in the upper 3 km below seafloor.To emphasise the differences, the percentage of change is calculated in Fig. 14c.
Close to the trench axis (CDP 25500-26000), a velocity reduction of more than 10 % from 2100 m/s to 1800 m/s is observed relative to the initial velocity (Fig. 14c).The final velocity close to this area increases from 1750 m/s at the seafloor to 2280 955 m/s at the plate boundary at 7400 m depth.In contrast, the uplifted sediment ridge (CDP 26000-27000) shows a velocity increase from 1750 m/s to 1850 m/s in the upper 500 m below the seafloor, whereas the velocity increases up to 2650 m/s at an observed basement high at 7100 m depth.If such a basement high hits the tip of the accretionary prism, the upper plate property would get fractured, deformed, uplifted at the leading edge of the subducting relief.While, as soon as the relief passes further down dip, the gravitationally driven subsidence and slumping of these sedimentary rocks will dominate at the trailing 960 edge of the basement high, and forms the escarpment at the basement high's trailing flank close to this trench location.This intensified structural modification of the upper plate by the subducting relief will result in a lossy and less consolidated rock status, compared to the major part of the accretionary wedge, and thus manifested in the local low-velocity zone inside the trench (Fig. 14b).
On the lower slope (CDP 27000-29000) an increase of the velocity of more than 10 % compared to the initial velocity is 965 observed and is comparable to the original OBS velocity.A thin pelagic layer of slope sediment with a maximum thickness of 100 m with velocities of 1750 m/s covers a not well structures accretionary prism.The velocity below the slope sediment increases gradually in the upper 1500 m up to 3400 m/s which is higher than the OBS velocity model.The relatively high velocity of the major part of the accretionary wedge, which is composed of the ancient oceanic pelagic sedimentary rocks, yields long-term compaction, consolidation of the sedimentary structure.Additionally, the complex reflectivity pattern of 970 strongly folded and fractured strata with limited spatial extents (Fig. 9a, 12d) due to compressional tectonic deformation manifests itself in small thrust ridges at the seafloor with landward dipping reflectivity pattern below (e.g.CDP 28300, CDP 28700).Deleted: Importantly, each velocity update to the velocity model will change the ray paths.The best solution is to slowly update the velocity in an iterative approach and not significantly disturb the pattern of the ray paths in the following iteration.Depending on the 1075 complexity of the geological structure and the approximated velocity deviation from the velocity model to the real subsurface geophysical system, the maximum allowed velocity update varies from case to case.In our cases, we set the maximum velocity update to 10-15 % of the current velocity model for each iteration, which is a compromise 1080 between computation time and result quality and stability.While a situation could occur where the local velocity update is more than 15 % of the current velocity model, a customized damping factor should be assigned to reduce the maximum velocity update, or the previous iteration needs to run twice.In comparison, if the update is smaller 1085 than 3 %, the inversion can be stopped or continued up to smaller grid size.A perfect final tomography result is never reachable due to the approximations of wave propagation.Qualitative control will give a comparison of the reflectors horizontal alignment in the CIP gather with respect to the previous iteration or initial iteration.To produce a 1090 final image, residual move-out correction of the smoothed NRM displacement field could always be applied directly to flatten all the remaining dipping events in the CIP gathers to maximize the quality of the migration CIP stacked image without additional changes to the velocity model (Perez and Marfurt, 2008).¶ 1095 4.2 Layer based vs. Grid-based Tomography ¶ Historically, two basic tomography schemes exist, the layer-based and the grid-based schemes, while hybrid schemes combine the advantages of both strategies.For a geological setting where a strong impedance interface limits the boundary of the medium property and 1100 the velocity change, a layer-based tomographic inversion with a prejudged bias information will often result in a good tomographic solution.The RMO of CIP-gathers is manually picked along strong continuous reflections defined as horizons to update the velocity for each layer in a top to bottom approach.For each layer, a vertical In between the dipping reflectivity patches, shallow deformed layered sediment structures with reduced velocities of 1900 m/s are observed with landward increasing thickness from 200 m to 500 m and starting to form anticline structures with the increasing spatial size and reflector continuity (CDP 28800-29000).
On the upper slope (CDP 29000-32400) the shallow reflector continuity from the lower slope increases in thickness from 1115 500 m up to 2000 m and forms continuous landward dipping structures (Fig. 9b, 10d).The folded anticline structures (CDP 29000-29600) at a depth of 4 to 6 km and a sequence of thrust ridges with intervening piggyback basin (CDP 30500 -31100, Moved up [7]: Once there is a conflict occurring between some of the equations in one grid cell, the minority picks (which could be good or bad) will be rejected by the tomography algorithm in order to 1120 get a stable and self-consistent result.Unrealistic picks have an unfavourable effect on the tomographic results when they become the majority.
Deleted: Several approaches should be applied to avoid this from happening, depending on the origin of the unreliable picks.For 1125 massive three-dimensional scattered crustal events below a basement structure, we digitized a horizon as a lower boundary of the available residual move-out picks (Fig. 4, 5a, and 5e).To attenuate side echoes, residual multiple energy, or cross-cutting events in the CIP-gathers (e.g., Fig. 2d, RMO_Pick_2) resulting in unexpected low velocities, 1130 we applied a dip-filter prior to the NRM displacement calculation.To separate regions of reliable and less reliable picks, the reflector coherency was calculated for the PSDM profile after each iteration.This coherency information was used as a weighting factor to the picks to distinguish whether the residual move-out originated from a 1135 well-stratified sediment environment or from a rough, discontinuous structure (Fig. 7d and 7f).¶

Grid-Based Tomography Result ¶
The Java accretionary margin is a localized end member of a foldand-thrust belt environment, with a complex geological background 1140 and small-scale structural heterogeneities.Along the lower and middle slope, strongly folded and fractured strata with limited spatial extents make the construction of a detailed layer-based macro model unachievable.Due to the high reflectivity in the upper 4 km below the seafloor, a grid-based tomography based on closely spaced (100 1145 m) CIP-gathers with pure data-driven depth error analyses provides the best solution.The initial subsurface velocity model (Fig. 11a) is a merge of manually heavily smoothed velocities estimated by interactive velocity picking and a smoothed version of the 2-D refraction seismic travel time tomography.The transition between the 1150 two velocity fields is marked by a semi-transparent grey line in Fig. 11a.Close to the trench axis and the uplifted sediment ridge (CDP 25500-27500), a velocity reduction of up to 10 % (Fig. 11c) indicates a low degree of compaction in the upper deformed sediment layers.Along the lower and upper slope (Fig. 11c, CDP 27500-32000), the 1155 opposite is observed.Aside from local sediment basins with a velocity decrease of up to 10 %, a significant velocity increase of more than 20 % in the upper 500 m indicates a longer history of compaction.Additionally, the complex reflectivity pattern (e.g., Fig. 9d) represents a long history of compressional tectonic deformation 1160 that manifests itself in thrust ridges at the seafloor and intervening sediment basins.The complex and laterally-disrupted reflectivity that has resulted from this history of deformation would make a conventional migration velocity analysis with sparse analysis locations extremely difficult.¶ To quantify the model uncertainties and to reduce mainly the migration computation time, new inversion strategies were 1240 developed by incorporating a Monte Carlo approach (e.g.Martin and Bell, 2019) and should be incorporated in the future.
Instead of getting one final model result, multiple model results were generated based on the sensitivity and resolution of the input data for the migration.To estimate the sensitivity and resolution, which is mainly determined by the acquisition parameters and the subsurface complexity, a checkerboard test with different wavelengths and magnitude of perturbation added to an initial velocity model will be inverted by a test tomography application.The difference to the initial model, namely the residual errors, can be used to constrain threshold values for model perturbations.The minimum spatial wavelengths and the maximum amplitude perturbation must be fulfilled for any velocity perturbation created randomly for a Monte Carlo simulation, but this analysis will not avoid the detection of side echo velocity anomalies.
To analyse model perturbations independently of the migration velocity, CIP-gather depth errors are de-migrated with their migration velocity.In the model domain, random perturbated populations of velocity input functions are generated, inverted and updated to the input velocity model to create a possible velocity model.All these populations of velocity models of the inversions can be statistically analysed, averaged and used for the next iteration of a migration.By analysing the cumulative depth error of the CIP-gathers after the migration iterations, convergence to a predefined obtained minimum depth error can be used to stop the inversion process automatically (Martin and Bell, 2019).

Anisotropic tomography
The wave propagation, namely direction and speed, is strongly influenced by the rock type and is generally depth and azimuth dependent.A fluid filled orientation of fractures or microcracks can cause anisotropy as well as a preferred orientation of minerals in the deeper crust.There exist several classes of symmetry for anisotropy.But for imaging and inversion often a simple transverse isotropic type with one axis of symmetry is assumed.The symmetry axis can be vertical (VTI), tilted (TTI) or horizontal (HTI).For a weakly anisotropic medium of acoustic waves, the dimensionless Thomsen parameters ε and δ are used to describe the ratio of the velocity variations (Thomsen 1986).
Complementary datasets like in this study of near-vertical reflections and wide-angle reflection and refractions with more horizontally propagating events offer one possibility to estimate the anisotropic parameters in the illuminated subsurface areas of both datasets.Classical modelling methods were nowadays replaced by inversion strategies due to the constant growth of observation density and increasing computational power.Several developments for weak anisotropy are published e.g.3D joint refraction and reflection tomography (Meléndez et al., 2019) or a ray-based gridded tomography for tilted TI media based on depth alignment of CIP gather (Wang and Tsvankin, 2011).
The ray-based gridded tomography (Eq.8 and 9) together with the non-hyperbolic NRM event tracking and picking can also be used to invert for the anisotropic parameters e.g.ε or δ.Based on an isotropic velocity and one Thomson parameter e.g.ε an initial anisotropic migration will be analysed in the CIP-domain and a depth error estimated.Instead of calculating a change in travel time, corresponding to a change of velocity / " (Eq.7), the calculation is modified to a change of the Thompson parameter /ε " .By exchanging  to  and solving Eq. ( 9), any CIP depth error is corrected due to a change of the parameter ε.An application to real data can be found in Woodward et al., 2008.The initial isotropic velocity should be ideally a velocity depth profile corresponding to a vertical seismic profile (VSP) at each location.To overcome this limitation a significant scale length smoothing  (Eq.9) needs to be applied as shown by (Wang and Tsvankin, 2011).
In the Java trench dataset, where near vertical and wide-angle and refracted OBS data both exist, a combined analysis is limited to the moveout sensitivity by the streamer length below 3-4 km of the seafloor.Additionally, both profiles do not completely coincide especially at the lower slope.The OBS model does not show significant velocity variations along the slope, especially not in the gravitationally driven slump area close to the trench axis.Instead, only a thin low-velocity layer with constant thickness is observed along the accretionary wedge (Fig. 6a).The data gap of OBS positions in the trench axis 1280 and the lack of incorporating local sediment basin into the initial wide-angle tomography model (Planert et al., 2010) have reduced the model reliability for further anisotropic analysis in the shallow illuminated area of both datasets.

Conclusions
The presented case study shows that CIP depth error estimations by depth warping in combination with a ray-based reflection tomography can improve depth-migrated images from MCS data.The non-rigid warping method provides reliable 1285 displacement fields for non-hyperbolic CIP depth errors.A semblance-based event tracking through the displacement field is limited by interfering reflected events.Due to the pure data-driven method of densely sampled depth error information (horizontal distance 100 m, vertical distance 50 m) more detailed spatial information for velocity corrections is available.In combination with a grid-based tomography, where depth errors are compensated by velocity changes, the inversion from longscale to short-scale lengths iteratively reduces the depth errors and improves the migration image.We suggest that further 1290 developments by integrating statistical analysis of the velocity updates (e.g.Monte Carlo approach) and extending the tomography for anisotropic parameters will provide new analyses tools for the subsurface image within the limits of ray-based methods.tomography with depth error estimations by Non-Rigid Matching (NRM) depth warping.Our results from analysis of NRM warping limitations in a synthetic CIP-gather, and a real data application from simple to complex sediment structures across the Java trench, revealed the following: ¶

1355
The NRM displacement field estimation in a relative referenced scheme between neighbouring traces is limited by coherent noise, intersecting events, or side echoes on 2D profiles.An alternative to NRM warping is a plane-wave destruction filter (PWD), but this method has reduced smoothness and accuracy.¶ 1360 The NRM method could be applied, as a purely data-driven method, to multi-channel seismic common image point gathers in the depth domain to estimate the residual depth error displacement field of any curvature.To calculate the relative depth error values in a gather along a predefined depth slice, a recursive depth variant correction 1365 followed by a cumulative summation of the individual displacements will yield the desired depth error information for the gather.Since the reflection depth tomography is capable of handling any kind of reflected arrival information without the restriction of curvature, the ... [10]

Figure 1 .
Figure 1.Processing scheme using Non-Rigid Matching in reflection tomography to update the velocity field during pre-stack depth migration of multi-channel seismic reflection data.The main processing steps are marked in red.

Deleted:Moved
An application of NRM for the vertical displacement calculation of neighbouring offset traces of a current 175 the current gather and the reference gather, which corresponds to the vertical smooth displacements between neighbouring traces of the current CIP gather 180 Deleted: gather depth error alignment Deleted: order curvatures (Sripanich must be depth variant shifted by the amplitudes of the second trace of the NRM field in Fig 2b to get aligned to the first trace.To further align the third trace of Fig. 2a, the 240 trace must be depth variant shifted by the amplitude of the third trace of the NRM field in Fig. 2b and additionally shifted the previous shifts which were applied to the second trace.The following equations are documented in the online repository (https://github.com/xyywy/Sup-SE2021-40)as web-based interactive script files.The depth correction alignment can be written as recursive formula Eq. (1): Figure 2. (a) Simulated complex geological 255 situations that would be frequently seen in pre-stack depth migrated (PSDM) common-image-point (CIP) domain.Deleted: A symmetrical diffraction, two interfering primaries with opposite polarity, and a non-linear local static undulation, including frequency versus offset signal variations.(b) The 260 Moved down [5]: NRM displacement of gather (a) calculated from trace n to the previous trace (n-1) for n > 1. (c) Application of the displacement correction from (b) to the gather of (a).Deleted: (d) Residual move-out picks calculated by recursive summation of the relative depth errors (b) at predefined depths to get 265 the cumulative depth error.¶ ¶ Deleted: of a trace to its previous trace Deleted: sections follows the general local dip trend well.An application of the NRM field to flatten the synthetic gather requires a 270 recursive depth variant correction.The NRM depth variant shifts of the first trace are applied to the first trace and all following traces.

Moved
Figure 2. (a) Simulated complex geological situations that would be frequently seen in pre-stack depth migrated (PSDM) common-imagepoint (CIP) domain.A symmetrical diffraction, two interfering parabolic events with opposite polarity, and parabolic event with a local curvature anomaly, including frequency versus offset signal variations.(b) The relative NRM displacement of gather (a) calculated from trace n to the previous trace (n-1) for n > 1. (c) Application of the displacement correction from (b) to the gather of (a).(d) Application of 335 [:]  represents the depth error relative to reflector depth  7 , array ℎ [:,:] represents the depth-corrected cumulative-summed NRM displacement field, and the index i represents the trace number.For any trace ℎ[",:]  in Fig.2dthe cumulative summed NRM field ℎ [:,:] in Fig.2fcan be calculated by: this situation are to either terminate permanently the current RMO-picks when Deleted: the crossing event or to terminate this picking procedure at the crossing point and start a new series of RMO-pick values again right after it, through to the end of the offset range.So, instead of 375 only one, two series or branches of depth error RMO picks Deleted: recorded and set as input information to the

Figure 3 .
Figure 3. (a) Comparison of residual depth error picks between the NRM and PWD method.(b) Depth difference between the depth errors of the NRM and PWD method.

Figure 4 .
Figure 4. (a) Simulated geological situations with multi sets of reflectors.(b) The NRM displacement of gather (a) calculated from trace n to the previous trace (n-1) for n > 1. (c) Application of the "relative-displacement correction" scheme from (b) to the gather of (a).(d) Residual move-out picks automatically calculated by the semblance-weighted grid-based approach of the relative depth errors (b).
contrast to the grid-based tomography, where vertical and horizontal velocity gradients are determined during the inversion, the layer-based tomography updates the lateral velocity variation between two user defined horizons with a predefined vertical velocity gradient.A comparison of layer-based and grid-based tomography results can be found in Riedel et al. (2019) andSugrue et al. (2004).Model areas, where velocities a priory are known, e.g. the water column above the seafloor, hybrid models are used to avoid tomographic velocity updates to propagate into selected areas.Furthermore, first order velocity contrasts, resulting in ray-path bending are problematic as they cannot be inverted by finite grid size.HereMoved (insertion) [7]hybrid models delivered the best results as shown byJones et al. (2007)  andFruehn et al. (2008).For most of the studies is common, that the grid-based tomography was applied to moderate layered structures in combination with the hyperbolic curvature scanning technique fromHardy (2003) to estimate the depth error in the CIP domain.An application result of a gridbased tomography combined with the NRM technique in a moderate layered structure is shown inCrutchley et al. (2020).3Application of a reflection tomography by CIP residual moveout warping across the Java trenchAn application of the NRM common image point depth error estimation in combination with an iterative grid-based tomography approach will be presented here in detail for three different structural settings along a profile crossing the Java trench.The complexity of the data examples increases from moderate horizontal layering, to dipping layered reflections, up to small disrupted dipping reflector elements.

Figure 5 .
Figure 5. Map of the study area offshore Southern Java.Local multi-beam bathymetry was acquired during SO190 cruise, overlain on the GEBCO_2020 grid (GEBCO, 2020).The location of the multi-channel and collocated wide-angle seismic profile is shown by a black line.Three examples at locations A, B, and C (marked in red) of the NRM based velocity updating for the depth tomography and pre-stack depth migration result are discussed in detail.Example A is located at shallow depth with simple complexity, whereas examples B and C are

515
Sequence Step NamesNormal and Nominal Geometry Establishment with CMP spacing of 6.25  m Anomalous and Random Noise Attenuation Padding Interpolated Traces to Zero Offset Interactive Velocity Analysis in Time Domain Initial Time-domain Velocity Building Shot Interpolation for Aliasing Elimination (from 50 m to 12.5 m shot distance) Surface-Related Multiple Prediction Multiple Attenuation 1: Frequency-Split 2D Cascaded Adaptive Filter Multiple Attenuation 2: Radon Dip Filter Multiple Attenuation 3: Inside Mute and Amplitude Clipping Kirchhoff Pre-Stack Time Migration Initial Depth Domain Velocity Building (Merge with Wide-Angle Model) * Kirchhoff Pre-Stack Depth Migration (PSDM) with Common Image Point Gather Output Pre-filtering of CIP (Common Image Point) Gather for NRM Calculation NRM Displacement Field Calculation CIP-Gather Residual Move-Out (RMO) Picks Calculation from NRM Field Dip Field, and Coherency Field Estimation from PSDM Section Map of the study area offshore Southern Java.520 Moved up [9]: The location of the multi-channel and collocated 540 wide-angle seismic profile is shown by a black line.Three examples at locations A, B, and C (marked in red) of the NRM based velocity Deleted: high complex structures where standard velocity analyses mostly fail by discontinuous highly dipping structures.

Deleted
collocated 2-D refraction seismic line covered by 46 OBS with a spacing of 6 km (Planert et al., 2010) and a velocity model 550

Figure 6 .
Figure 6.(a) The original OBS velocity model with the line drawing based on the final PSDM image.(b) The initial velocity model for the reflection tomography merged from the multi-channel seismic velocity analysis above the white transparent band and the wide-angle velocity model (below the white transparent band).The line drawing is based on the final PSDM image.580 To balance the quality of the picks, Moved down [10]: ranging from CDP 46700 to CDP 50800 is displayed in Fig.4Deleted: 4a.The resulting initial Kirchhoff pre-stack depth migration (Fig.4c660 Moved down [11]: ) retrieves a coherent image of the shallow sedimentary portion, while the energy in the deeper part close to the basement is not very well collapsed, resulting in a series of overmigrated events.The displayed reflector dip field (Fig.Deleted: 4e) and coherency field (Fig.4f) are extracted from the 665 final migration section.

Figure 7 .Figure 8 .
Figure 7. Depth tomography example A from Fig. 5, with CDP ranging from 46700 to 50800.(a) Initial velocity model merged from velocity analysis and wide-angle refraction tomography.(b) Final velocity model after five iterations of NRM based depth tomography and PSDM.(c) PSDM result based on the initial velocity model.(d) PSDM result based on the final velocity model.(e) Reflector dip field calculated 670

Figure 6 .
Figure 6.Depth migration stack section of six Moved down [14]: iterations depth tomography with final velocity model overlain.The location of the profile is illustrated in Fig. 3 Deleted: 3

Figure 9 .
Figure 9. Depth migration stack section of five iterations depth tomography with final velocity model overlain.The location of the profile is illustrated in Fig. 5 as a yellow line.Rectangular boxes are discussed in the section on field data examples from the accretionary wedge.

Figure 10 .
Figure 10.Depth tomography example B in Fig. 5 and Fig. 9b, with CDP ranging from 29300 to 30500.(a) Initial velocity model merged from velocity analysis and wide-angle refraction tomography.(b) Final velocity model after five iterations of NRM based depth tomography and PSDM.(c) PSDM result based on the initial velocity model.(d) PSDM result based on the final velocity model.(e) Reflector dip field calculated from the final PSDM result.(f) Reflector coherency field calculated from the final PSDM result.Notice the continuity and reflector dip change of the folded sediment layers at a depth of 4.0-4.8km in (c) and (d) based on the change of the initial velocity and final velocity 825

Figure 11 .
Figure 11.NRM velocity updating of Fig. 10 in CIP domain.(a) CIP gathers based on the initial velocity model.(b) CIP gathers based on the final velocity model.(c) Initial NRM depth shifts in the CIP domain.(d) Final NRM depth shift in CIP domain.(e) RMO picks calculated from the initial NRM displacement field.(f) RMO picks calculated from the final NRM displacement field.Notice that the distinct area of 830

Figure 12 .Figure 13 .
Figure 12.Depth tomography example C in Fig. 5 and Fig. 9a, with CDP ranging from 27200 to 28600.(a) Initial velocity model merged from velocity analysis and wide-angle refraction tomography.(b) Final velocity model after six iterations of NRM based depth tomography 930 depth errors distributed over the entire offset range yields the estimation of the interval velocity changes along their ray paths between the source Deleted: receiver (Jones, 2010).In grid-based tomography, the updated subsurface velocity model is smoothed by a predefined 985 scale length that decreases in size during the iterations (Tab.2).To ensure a stable result over all scale lengths, each iteration started with the largest scale length application sequence one up to the desired target scale length iteration sequence.Because the starting velocity was estimated based on a wide-angle 990 tomography (Planert et al., 2010), the first iteration with a largescale length did not experience significant velocity corrections.The higher the iteration number with corresponding reduced scale length, the more pronounced were the velocity updates observed.After the sixth iteration, the velocity updates were 995 stopped because in subsequent iterations, the fluctuations of velocity changes could no longer be related to subsurface structures.

1105
Figure 14.(a) The initial velocity model merged with the wide-angle velocity model (below the white transparent band).The line drawing is based on the final PSDM image.(b) The final velocity model calculated from five iterations of the ray-based tomographic velocity 1110

Figure 15 .
Figure 15.(a-e) Velocity update after each iteration of the gridded tomography together with the CIP depth alignment after each iteration.(f) The cumulative sum of the velocity updates from all five iterations together with the horizontal and vertical smoothing length used by the scale length reductions in each iteration.
We have successfully performed a grid-based reflection 1350 Based on the above discussion, any intermediate irregular index  [ 2 ] in the array  [ : ] can be simply expressed by  [ 2 ] =  − ℎ [ ",-] , where the j represents the original regular depth index of the array  [ ",: ] , and the k represents the index of the NRM

Table 1 .
Seismic processing sequences and image grid sizes.

Table 2 .
Regional depth variant water velocity extracted from the Climatological Atlas of the World Ocean MB-System.545

Table 3 .
Successive smoothing scale length reduction applied for each iteration of the depth tomography.