Geophys. J. Int. (2009) 179, 201–214

doi: 10.1111/j.1365-246X.2009.04213.x

GJI Seismology

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

2-D/3-D multiply transmitted, converted and reflected arrivals in complex layered media with the modified shortest path method
Chao-Ying Bai,1,2 Xiao-Ping Tang1 and Rui Zhao1
1Department of Geophysics, College of the Geology Engineering and Geomatics, Chang’an University, Xi’an 710054, China. E-mail: chaoying_bai@yahoo.com 2Key Laboratory of Western China’s Mineral Resources and Geology Engineering, Education Ministry of China, Xi’an 710054, China

Accepted 2009 April 14. Received 2009 April 12; in original form 2008 December 18

SUMMARY Grid-cell based schemes for tracing seismic arrivals, such as the finite difference eikonal equation solver or the shortest path method (SPM), are conventionally confined to locating first arrivals only. However, later arrivals are numerous and sometimes of greater amplitude than the first arrivals, making them valuable information, with the potential to be used for precise earthquake location, high-resolution seismic tomography, real-time automatic onset picking and identification of multiple events on seismic exploration data. The purpose of this study is to introduce a modified SPM (MSPM) for tracking multiple arrivals comprising any kind of combination of transmissions, conversions and reflections in complex 2-D/3-D layered media. A practical approach known as the multistage scheme is incorporated into the MSPM to propagate seismic wave fronts from one interface (or subsurface structure for 3-D application) to the next. By treating each layer that the wave front enters as an independent computational domain, one obtains a transmitted and/or converted branch of later arrivals by reinitializing it in the adjacent layer, and a reflected and/or converted branch of later arrivals by reinitializing it in the incident layer. A simple local grid refinement scheme at the layer interface is used to maintain the same accuracy as in the one-stage MSPM application in tracing first arrivals. Benchmark tests against the multistage fast marching method are undertaken to assess the solution accuracy and the computational efficiency. Several examples are presented that demonstrate the viability of the multistage MSPM in highly complex layered media. Even in the presence of velocity variations, such as the Marmousi model, or interfaces exhibiting a relatively high curvature, later arrivals composed of any combination of the transmitted, converted and reflected events are tracked accurately. This is because the multistage MSPM retains the desirable properties of a single-stage MSPM: high computational efficiency and a high accuracy compared with the multistage FMM scheme.
Key words: Numerical solutions; Seismic tomography; Computational seismology; Wave propagation.

1 INTRODUCTION
Multiply transmitted, converted and reflected (or refracted as a special case) seismic phases are caused by both continuous and discontinuous variations in wavespeed within the crust and the Earth’s deeper interior (e.g. the Moho discontinuity). These complex arrivals are recorded by modern digital broad-band seismographs worldwide, and they include more subsurface information than the first (direct) arrivals alone. However, traditionally in earthquake location, traveltime tomography and other applications, only the information of first arrivals (or other specific seismic phases) is used. Nonetheless, later arrivals are important for improving the location accuracy (e.g. Kennett & Engdahl 1991) and also tomographic resolution (e.g. Zhao et al. 2005) because more constraint parameters are involved in seeking a physical solution. The problems of ex-
C 2009 The Authors Journal compilation C 2009 RAS

ploiting information from complex arrivals on digital broad-band seismograms are that (1) there is a lack of robust algorithms to trace such multiple (or later) seismic phases with an appropriate computational efficiency and solution accuracy and (2) automatic multiple seismic phase detection and recognition on digital broad-band seismograms is still under development (e.g. Di Stefano et al. 2006; Reading et al. 2008). In this study, we focus on how to generate a variety of multiple seismic phases within a complex layered velocity model, when a relatively high degree of interface curvature is involved. Our study provides a step toward being able to exploit fully the information obtained by the digital broad-band seismographs using a practical algorithm.
There are many techniques available for tracing (high frequency) seismic arrivals in laterally heterogeneous media, but most are only capable of yielding the first arrivals and/or single (or primary)
201

202 C.-Y. Bai, X.-P. Tang and R. Zhao
reflected phases between source and receiver. Traditionally, the ‘shooting’ and ‘bending’ methods are used, as described by Julian & Gubbins (1977), Um & Thurber (1987) and Ce˘rveny (1987, 2001). In the shooting method, ray tracing is solved as an iterative initial value problem. Shooting methods can be very efficient at tracing seismic rays if there is a large set of rays and closely spaced receivers. The major shortcoming of the method is that the target points are frequently not well-behaved functions of the shooting angle, making some rays impossible to trace (e.g. receivers located in the shadow zone). In the bending method, ray tracing is handled as an iterative boundary-value problem. Bending methods are advantageous over shooting methods when the receiver positions are ill-behaved functions of the shooting angle. The drawbacks of both of the above ‘two-point’ ray tracing algorithms are as follows: (1) the solutions sometimes can be trapped in a local minimum and (2) it is necessary to repeat the ray tracing when the location of the receiver is changed. This is, therefore, extremely costly in CPU time, making such an approach impractical, for example, for large 3-D models.
Other approaches for seismic ray tracing have been devised, such as the Lagrangian solution to the problem of seismic ray tracking, introduced by Lambare´ et al. (1992) and Vinje et al. (1993). Here a set of points is used to represent the seismic wave front, which is iteratively evolved by using local ray tracing with a given time step. Other methods, originally developed outside of geophysics, include the essentially non-oscillatory schemes (e.g. Harten et al. 1987; Shu & Osher 1989), and the weighted essentially non-oscillatory schemes (e.g. Liu et al. 1994; Jiang & Shu 1996; Jiang & Peng 2000). These methods exploit algorithms for solving Hamilton– Jacobi type equations. The ‘phase space’ method (Osher et al. 2002; Fomel & Sethian 2002; Coman & Gajewski 2005) can also be used for tracking multivalued seismic ray paths (e.g. Hauser et al. 2008).
In the later part of the 20th century, grid-cell based algorithms became popular because they offer a global solution with robust simulation and yield more accurate solution and consume less CPU time. Two practical approaches, in particular, have shown promising results for first-arrival calculations. The first method is the finite difference (FD) eikonal equation solver (Vidale 1988, 1990, and revisions thereafter, e.g. Van Trier & Symes 1991; Cao & Greenhalgh 1994; Schneider 1995; Afnimar & Koketsu 2000; Kim 2002; Qian & Symes 2002). Several studies have concentrated on how to use eikonal equation solvers to compute reflected (or refracted) traveltimes in layered media (e.g. Podvin & Lecomte 1991; Riahi & Juhlin 1994; Hole & Zelt 1995). With the advent of the fast marching method (FMM, Sethian 1996), Rawlinson & Sambridge (2004) developed a multistage scheme to trace multiply reflected and transmitted arrivals in complex 2-D layered media with the FD eikonal equation solution. Later, this multistage scheme was extended to three dimensions to trace multiple reflections (or refractions) in heterogeneous media (e.g. De Kool et al. 2006). Such improvements significantly elevate the status of grid-cell based algorithms in seismic ray tracing, where multiple arrival information has been urgently needed.
The other method is referred to as the shortest path method (SPM) (e.g. Nakanishi & Yamaguchi 1986; Moser 1991; Cao & Greenhalgh 1993; Fischer & Lees 1993; Klimesˇ & Kvasnicˇka 1994; Cheng & House 1996; Gruber & Greenhalgh 1998; Zhang & Toksoz 1998; Van Avendonk et al. 2001; Bai et al. 2007). To avoid solving a differential equation, a network (or graph) is formed by connecting all neighbouring nodes to offer many possible seismic ray paths. The key issue here is then to find the ray paths with the minimum traveltimes (Fermat’s principle) between the source and all nodes within

the network (similar to a road map where we find the most direct route between two places). The modified SPM (MSPM) approach (Bai et al. 2007) allows us to trace complex arrivals in layered media. It is based on the idea of a multistage scheme (Rawlinson & Sambridge 2004) for tracing multiply reflected, transmitted and converted seismic phases. Each time a wave front interacts with an interface, the onward propagation of that wave front is treated as a new computational domain.
In this paper, we conduct a comparison between the multistage FMM and the multistage MSPM in terms of solution accuracy and computational efficiency. We then present several examples to assess the stability and performance of MSPM in simple and complex 2-D/3-D layered media (including the Marmousi model). We aim to demonstrate that the multistage MSPM approach retains the desirable characteristics of the single-stage MSPM and is a practical way of tracing multiple arrivals in layered media.
2 METHODOLOGY
First, we will briefly review some approaches used in tracing the primary (or single) reflection and then outline the multistage scheme implemented with the MSPM approach.
2.1 Two-way approach
To find primary reflections, we may use either the FD eikonal equation solver or regular SPM approach, where the first-arrival traveltime fields from both the source and the receiver to the entire sampled interface are calculated (hereafter referred to as the twoway approach). Then Fermat’s principle of stationary time is used to locate the reflector along the interface. One problem with the twoway approach is that two or more reflectors may be located when a strongly curved (or irregular) interface is involved. To guarantee finding a single reflector for a specified source–receiver pair, one should introduce Snell’s law to constrain the reflector.
For 3-D applications there are two ways to proceed: (1) assume that there are no scatterers in the ray propagation field, which means that the reflectors are constrained to the source–receiver plane, and in such a case, the 3-D problem reduces to the 2-D case for locating the reflector. (2) If the above assumption is invalid, it means that the ray is capable of being scattered, and the Fermat’s principle of stationary time together with the Snell’s law should be used at all sampled nodes in the subsurface interface. This increases the computational costs in proportion to the number of sampled node spacing on the subsurface.
Using the two-way approach, a traveltime field needs to be computed for each source and receiver, which adds an extra amount of CPU time. Furthermore, it is difficult to simulate multiple times of reflected, converted and transmitted phases between the interfaces. For these reasons, we have to seek an alternative approach.
2.2 The one-stage MSPM
The MSPM differs from the regular SPM in that (1) primary and secondary nodes are introduced in the model parametrization (refer to Fig. 1), (2) a bilinear (2-D case) or trilinear (3-D case) velocity interpolation function is used to calculate the velocity values at the secondary node positions (including the source and receiver locations if they do not coincide with the primary node positions), (3) there is no need to use the forward star technique as the regular SPM
C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Multiple arrival tracking with MSPM 203

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Figure 1. A diagrammatic illustration for computing the transmitted, converted and reflected arrivals in multistage MSPM approach.

does. The forward star technique refers to the geometric arrangement of network connections or possible ray branching points into adjacent nodes. (for detailed explanation, see Klimesˇ & Kvasnicˇka 1994), because the outward extended wave front computation is accomplished cell by cell.
In calculating the minimum traveltimes and locating the related ray paths for all the grid nodes, we can gradually expand the volume of the computed nodes by continually adding the undetermined neighbouring nodes to the computed nodes. In this process, one should start with the node that has a minimum traveltime in the subset N j (where N j is the total number of computed nodes in the current wave front) to keep track of the first arrival times for the undetermined nodes. An interval sorting method (Klimesˇ & Kvasnicˇka 1994) is used so that the longer traveltimes are deleted; only the minimum traveltime and the related ray path are recorded. The minimum traveltime from a source node i to an undetermined node j in a cell is expressed as

tij

=

min
i∈N j

ti

+

2D(xi , x j ) [V (xi ) + V (x j )]

,

(1)

where D(xi , x j ) is the distance between the source node i and the undetermined node j and V (xi ) and V (x j ) are the velocity values at the ith and the jth node positions. If the ith or the jth node does

not belong to the subset of primary nodes, then the velocity value

at such a node position in a cell can be obtained by a Lagrangian

interpolation function

⎛

⎞

V (x)

=

8 k=1

⎜⎝ 8
l =1

(x (xle

− xek ) − xek )

⎟⎠V

(xek

),

(2)

l =k

where xke and V (xke) (k = 1, 2, ......, 8) are the vector coordinates and the sampled velocity values at the primary nodes of the cell, respectively. Meanwhile, in this process, the order number of the incident point (node) i∗giving the minimum traveltime to the node j is recorded for the coordinates of the related ray path.
The accuracy of the MSPM is comparable to the accuracy obtained from the regular SPM. Furthermore, MSPM has the follow-

ing advantages over the regular SPM (for more detail, see Bai et al. 2007):
(1) Less computer memory and CPU time are used (typically one order of magnitude).
(2) Easy handling of media containing large velocity contrasts (Bai & Greenhalgh 2005a). This is due to velocity values being bilinearly or trilinearly linked to secondary nodes or source and receiver positions.
(3) Fewer number of total cells per simulation (roughly equal to the primary nodes), which leads to easy incorporation with any inversion subroutines to solve real tomographic problems (Bai & Greenhalgh 2005b; Bai et al. 2008).
2.3 The multistage MSPM
Figure 1 gives a schematic representation of the multistage MSPM scheme. After dividing the model into several layers corresponding to the velocity structure (for example, the upper and lower computational domains separated by the interface shown in Fig. 1), a simulated downwind wave front (P1 phase) is propagated through the upper computational domain (or active region) until it impinges on all sampled nodes of the interface (In our notation convention: P or S represents a P or S wave, respectively, and the number of subscript or superscript indicates a downwind or upwind seismic wave propagating in different regions, respectively.). At this stage the independent computational domain is halted at the active region and we are left with a narrow band of traveltime values defined along the sampled interface. From here, a downwind propagation of a transmitted wave (P1 P2 phase) or transmitted and converted branch (P1 S2 phase) can be simulated by reinitializing it, starting at the sampled node position with the minimum traveltime (i.e. Huygens principle, the node is treated as one new source point in the wave front). The wave is then projected into the adjacent (lower) layer (the waiting region). Meanwhile, an upwind-propagating wave front consisting of a reflected branch (P1 P1 phase) or a reflected and converted branch (P1 S1 phase) can now be obtained by reinitializing wave front and starting at the sampled node position with

C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

204 C.-Y. Bai, X.-P. Tang and R. Zhao

2006). Because the code of the multistage FMM is in Cartesian coordinates for 2-D models and spherical coordinates for 3-D models, but our code for multistage MSPM is in Cartesian coordinates under both 2-D and 3-D models, the following comparisons are based on the 2-D velocity model. Additionally, the multistage FMM is only capable of tracking the multiply transmitted and reflected arrivals, and therefore, we study multiple arrivals of transmissions and reflections.

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Figure 2. The traveltime errors against the analytic solutions along the receiver line for the multistage FMM solutions (diagram a) and the multistage MSPM solutions (diagram b) with changing grid spacing as indicated in the figure.
the minimum traveltime, from the narrow band of the interface into the incident layer. Different velocity models (i.e. P or S) are used if wave mode conversion occurs at the interface. Fig. 1 illustrates the four basic branches of seismic wave propagation off the velocity interface. Multiple arrivals composed of any combination of these four events can be built up by using this method. In summary the multiple arrivals are the different combinations or conjugations, via velocity discontinuities (i.e. interfaces), of the incident, reflected and converted phases, which obey Snell’s law, Fermat and Huygens principle.
3 BENCHMARK TESTS WITH THE FFM SCHEME
Previously we noted that for first arrival calculations, the one-stage MSPM delivers the same solution accuracy under less CPU time than the regular SPM approach (generally it is one order of magnitude faster—see Bai et al. 2007). To assess the viability of the multistage MSPM, we compare it against the multistage FMM (2D case, Rawlinson & Sambridge 2004; 3-D case, De Kool et al.

3.1 Primary reflections in a uniform velocity field
First, we compare the algorithms when the same grid spacing is employed in the model parametrization. We select a constant velocity model (4.0 km s–1) with a scale length of 100 km × 40 km. In the model parametrization we select four different grid spacing (that is 1000, 500, 250 and 125 m) for the FMM approach. For the MSPM approach, 3, 7, 15 and 31 secondary nodes are added to each cell boundary (cell size is 4.0 km × 4.0 km) to keep the corresponding node spacing (note that there are no nodes at all inside the cell for the MSPM). In the multistage FMM approach the second order and source refined scheme was used as Dr Nicholas Rawlinson suggested (Rawlinson & Sambridge 2004; Rawlinson, personal communication). The source was placed at the upper lefthand corner of the model and 100 receivers with 1.0 km uniform spacing were located along the top boundary of the model, and a horizontal reflecting line was located at the depth of 30.0 km. The traveltimes of the primary reflection were calculated at each receiver location for both the FMM and the MSPM approaches. Fig. 2 shows the traveltime errors against the analytic solutions along the receiver line for the FMM approach (diagram a) and the MSPM approach (diagram b). The average time and relative errors, and the corresponding CPU times are given in Table 1 (the PC computer we used here is Dell Optiplex 755–2.53 GHz). From Fig. 2 and Table 1, we see that except for the case of three secondary nodes added in each cell boundary (such selection is never used in the MSPM approach, see Bai et al. 2007), the MSPM approach is more accurate than that of the FMM in computing the primary reflections. For the MSPM algorithm the ray paths of node connections (numerical solutions) coincide exactly with or very close to the ray path of the analytic solution in some directions, so that the traveltime errors in such directions are zero or very small. Hence the errors along the receiver line fluctuate within a specific error bar. Meanwhile in the FMM algorithm the errors increase slightly with the source–receiver distances but are fairly flat in each case. Such features will be discussed further in later tests. From Table 1, we see that the CPU times are several times faster in the MSPM approach compared with the FMM approach if the same grid spacing is used to parametrize the velocity model.

Table 1. Summary of the averaged time and relative errors against the analytic solutions for both the multistage FMM and the MSPM approaches and the related CPU time consumption in tracking the primary reflections in uniform velocity field with changing grid spacing.

Grid spacing (m)

Total number of nodes

FMM

MSPM

Time error (ms)

FMM

MSPM

Relative error (per cent)

FMM

MSPM

CPU time (s)

FMM

MSPM

1000 500 250 125

4141

2073

53.9

76.7

0.27

0.39

0.3

0.2

16 281

4313

26.0

20.0

0.13

0.10

0.7

0.3

64 561

8793

12.8

5.7

0.06

0.03

2.4

0.7

25 7121

17 753

6.3

1.4

0.03

0.01

8.6

1.9

C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

Multiple arrival tracking with MSPM 205

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Figure 3. Two 1-D linear velocity models (indicates the corresponding ray paths) used in the benchmark tests for investigating the solution accuracy and computational efficiency (diagram a, Model-1 and diagram b, Model-2). In figure the different grey colour regions indicate the different computational domains, the white horizontal lines are the reflected lines and the black lines are ray paths.
3.2 Multiple transmissions and reflections in a 1-D linear velocity field
In this test, two 1-D linear velocity models (with vertical velocity gradient 0.04 km s–1) were selected (see Fig. 3). In Model-1 (Fig. 3a), two horizontal reflectors are located at the depths of 15 and 30 km, and in Model-2 three horizontal reflected lines are located at the depths of 10, 20 and 30 km, respectively. In both models the source is again placed at the top left-hand corner of the model, and 100 receivers are located along the top boundary of the model with 1.0 km uniform spacing. The traveltimes for transmitted and reflected arrivals were computed at each receiver location (for the corresponding ray paths, see Fig. 3a for Model-1 and Fig. 3b for Model-2). Fig. 4 shows the time difference between the FMM and the MSPM approaches [i.e.FMM(t) − MSPM(t)] along the receiver line (diagram a, results of Model-1 and diagram b, results of Model-2). The averaged time differences, the averaged relative values [i.e.|FMM(t) − MSPM(t)| × 100.0/FMM(t)] and the corresponding CPU times are listed in Table 2 (results of Model-1) and Table 3 (results of Model-2). From Fig. 4 and Tables 2 and 3, we see that the traveltime differences between the FMM and the MSPM are increased slightly with increasing source–receiver distance. This is partially due to the behaviour of the FMM approach, in which the traveltime errors increase slightly with the source–receiver distance (see Fig. 2a). However, the other cause is that the MSPM is capable of improving the accuracy of the computed traveltimes much faster than that of the FMM in the grid spacing reducing process. To reinforce this point, we present the T–D curves for both the FMM and the MSPM approaches for computing the Model-2 using the finest grid spacing (Fig. 5). From Fig. 5 it is evident that the numerical solutions of the MSPM are systematically smaller than that of the FMM, and the T–D curve of the MSPM is much smoother than that of the FMM. As we already know, the numerical solutions of both the FMM and the MSPM approaches are larger than or sometimes equal to the analytic solutions in the uniform or 1-D linear velocity
C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

Figure 4. The traveltime differences along the receiver line between the FMM and the MSPM approaches in computing the Model-1 (diagram a) and the Model-2 (diagram b) with indicated varied grid spacing.
field. In this respect, the smaller the numerical solution is, the more accurate the algorithm.
In summary, for the multistage scheme to track the multiple arrivals in 2-D layered media, both the solution accuracy and the computational efficiency of the MSPM approach are superior to the FMM approach if the same grid spacing is employed in the model parametrization. Such features can be reasonably expected also in the computation of the 3-D model (see Bai et al. 2007).
4 SIMULATION EXAMPLES
After conducting the above comparisons for both the solution accuracy and the computational efficiency, now we test a relatively simple 2-D layered velocity model to see how the multistage MSPM scheme can be used to simulate multiple arrivals of the transmitted primary reflections. Next, we investigate the validity of the multistage MSPM against the Marmousi velocity model. A number of complicated combinations of the transmitted, reflected and converted phases are then located, and finally a 3-D application is performed to show the robust and global solution of the multistage MSPM approach.
4.1 Transmitted primary reflections in the 2-D layered model
We now consider only the transmitted primary P-wave reflections. The P-wave model can be divided into three independent computational domains (referred as layers 1, 2 and 3 from top to bottom) with three curvilinear interfaces (named as interfaces 1, 2 and 3) separating them (see Fig. 6). The velocity fields are specified in the figure, and we effectively have an undulating-layered model. For clarity in Fig. 6, only one secondary node (M = 1) is added onto each cell boundary (actually more secondary nodes are needed for high accuracy, see Bai et al. 2007). In the model parametrization, the

206 C.-Y. Bai, X.-P. Tang and R. Zhao

Table 2. Summary of the averaged time differences and the relative values between the multistage FMM and the MSPM approaches, and the corresponding CPU time consumption in tracking the transmitted and multiple reflections (ray paths see Fig. 3a) for in 1-D linear velocity field with varied grid spacing.

Grid spacing x FMM

Nodes of FMM

Nodes of MSPM

Time difference (ms)

Relative value (per cent)

CPU time of FMM (s)

CPU time of MSPM (s)

1000 500 250 125

4141

2073

40.8

16 281

4313

54.0

64 561

8793

63.1

257 121

17 753

63.6

0.16

0.4

0.2

0.21

1.0

0.5

0.24

3.3

1.0

0.25

12.0

2.8

Table 3. Summary of the averaged time differences and the relative values between the multistage FMM and the MSPM approaches, and the corresponding CPU time consumption in tracking the multiple transmissions and reflections (ray paths see Fig. 3b) in 1-D linear velocity field with changing grid spacing.

Grid spacing (m)

Nodes of FMM

Nodes of MSPM

Time difference (ms)

Relative value (per cent)

CPU time of FMM (s)

CPU time of MSPM (s)

1000 500 250 125

4141 16 281 64 561 257 121

2483 5163 10 523 21 243

58.8 124.9 146.1 150.0

0.16

0.5

0.4

0.34

1.5

0.8

040

4.7

1.9

0.41

17.1

5.1

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

interfaces 2 and 1 are depicted in Figs 7(d–f). Fig. 8 shows the entire wave front propagation starting from the source and transmitting through interfaces 1 and 2, and reflecting back at interface 3 and transmitting again through interfaces 2, 1 and towards the top surface of the model.

Figure 5. The T–D curves in computing the Model-2 for both the FMM and the MSPM approaches with the finest grid spacing.
length of the rectangular cell is 1.0 km and the number of secondary nodes (M) added onto each cell boundary is 9. Therefore, the node spacing along the cell boundaries is 100 m. The sampled node spacing on the interface is also the same as that of the cell boundaries. The source was set at the upper left-hand corner of the model and the downwind-propagating wave fronts transmitted through interfaces 1 and 2 are displayed in Figs 7(a–c). The upwind-propagating wave fronts reflected back from interface 3 and transmitted again through

4.2 Multiply transmitted and reflected arrivals in the Marmousi model
It is important to test the new procedure in a more complex velocity field. For this purpose the Marmousi model was selected to see how the multiple seismic waves propagate in a heterogeneous medium. Two interfaces were introduced according to the velocity features of the hard Marmousi model (suppose that the bottom boundary of the model is another horizontal interface, see Fig. 9). In the model parametrization the length of the rectangular cell is 24 m and the number of secondary nodes (M) added onto each cell boundary is 15. This time the source was set at the centre of the top surface of the model. Again Figs 10(a–c) displays the simulated wave fronts emanating from the source, propagating downwind and transmitting through interfaces 1, 2 and to the bottom interface. The upwind wave fronts reflected back from the bottom horizontal interface and transmitted through interfaces 2 and 1 are captured in Figs 10(d–f). Fig. 11 is the entire wave front propagation starting from the source and transmitting through interfaces 1, 2 and toward the bottom boundary (Fig. 11a), and reflecting back from the bottom interface and transmitting again through interfaces 2, 1 and toward the top surface (Fig. 11b). Comparing the wave front propagation between the layered model and the Marmousi model, more complex wave fronts should be expected in the latter, as in Fig. 11. The wave fronts are enlarged along the high-velocity zones and flattened along the low-velocity zones. Meanwhile, the wave fronts deflect when they pass through a velocity discontinuity (see Figs 11a and b).
4.3 Reflected, transmitted and converted phases in the 2-D layered model
The above examples only deal with reflections and transmissions of the same wave types (e.g. P wave), but the real situation is more
C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

Multiple arrival tracking with MSPM 207
Figure 6. Parametrized 2-D undulating-layered P-wave velocity model with three curvilinear interfaces (black curve lines). In the figure: black circles, primary nodes; grey circles, secondary nodes and white circle, samples nodes on the interfaces. Note that the model can be divided into three independent computational domains (or layers), separated by the interfaces.

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Figure 7. The downwind wave fronts propagating through the first layer (a), secondary layer (b), the third layer (c) and upwind wave front reflected back from the interface 3 and propagating through again the third layer (d), the secondary layer (e) and the first layer (f) of the Model of Fig. 6. The source was set at the top left-hand corner and the wave front interval is 0.4 s.

Figure 8. The entire downwind wave front propagation (dark curves) and upwind wave fronts (thin curves) reflected back at the third interface. The wave front interval is 0.5 s. In the figure, the black circle is the source and dark curves are the interfaces.
C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

complex. There also exist mode conversions from one wave type to another (e.g. P to S or S to P) when the seismic rays penetrate through or echo back from the velocity discontinuity (or interface). Here, we only discuss the multiple transmitted and converted primary reflections. We still use the same model of Fig. 6, where the source position was set at the top left-hand side (X = 5.0 km) and 11 receivers are placed along the top right-hand side (starting from X = 45.0 to 95.0 km with 5.0 km spacing interval). We suppose that the S-wave model has the same structures and interfaces as the P-wave model, but with reduced velocity value. That is the S-wave model is obtained from the P-wave model by dividing a factor of 1.732 (VP / VS = 1.732).
Figure 12 shows the ray paths of multiply transmitted and converted primary reflected phases P1 S2 P2 P1 (Fig. 12a) and P1 P2 S2 P1

208 C.-Y. Bai, X.-P. Tang and R. Zhao
Figure 9. Marmousi model with the two interfaces (black curves) and basal interface (black line). Note that the two curved interfaces were drawn according to the velocity variations.

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Figure 10. The downwind wave fronts propagating through the first layer (top left-hand), secondary layer (middle left-hand), the third layer (bottom left-hand) and upwind wave front reflected back from the model bottom boundary and propagating again through the third layer (bottom right-hand), the secondary layer (middle right-hand) and the first layer (top right-hand) of the Marmousi model. The source was set at the top surface centre and the wave front interval is 0.05 s. The figures are superposed onto the Marmousi model.

(Fig. 12b). Obviously the P1 S2 P2 P1 and P1 P2 S2 P1 phases transmit and convert when crossing the interface 1, or convert when reflecting at the interface 2 in a different way. Fig. 13 displays the ray paths of multiply transmitted and converted primary reflected phases at interface 3 (diagram a, P1 S2 P3 S3 P2 S1 phase and diagram b, P1 P2 S3 P3 S2 P1 phase). The seismic rays propagate very differently when transmitting through, converting and reflecting at the different interfaces, dependent on the velocity variation and interface undulations.
4.4 Multiple transmissions, reflections and conversions in the 3-D layered model
We find that there is no significant difference for the multistage MSPM in 2-D or 3-D applications, but it will take more CPU time with a 3-D model, because many cells and nodes will be

created in the 3-D model parametrization. Theoretically, one of the advantages of the multistage MSPM when applied to the 3-D case is that it is capable of simulating seismic scattering. The 3-D velocity model was previously used by Bai and Greenhalgh (2005b) to test a tomographic inversion procedure to recover the fine 3-D velocity structure of the crust. Our model, shown in Fig. 14, has horizontal dimensions of 40 km × 40 km and a vertical extent of 30 km. The velocity distribution is smooth, and fixed at 3 km s–1 everywhere along the top surface of the model. At the base of the model (30 km depth) the velocity oscillates between 6.15 and 7.35 km s–1, with the maxima and minima 10 km apart. The velocity value in a box set within the horizontal plane either increases (low-velocity zones) or decreases (high-velocity zones) from its centre outward (see Fig. 14a for details). The velocity in the vertical direction increases linearly with depth (see Fig. 14b for details). The background velocity of the model base (Z = –30.0 km) is 6.75 km s–1. The variable velocity patterns are the same at the different depths, but the scale of the

C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

Multiple arrival tracking with MSPM 209
Figure 11. The entire downwind wave fronts propagating through the Marmousi model via the first, second and toward the bottom interfaces (upper diagram) and the wave fronts, reflected back from the bottom interface, propagating upwind (via interfaces 2 and 1) toward the top surface (lower diagram). Both figures are superposed onto the Marmousi model and the wave front interval is 0.075 s.

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Figure 12. Ray paths of transmitted and converted primary reflected arrivals P1S2P2P1 (diagram a) and P1P2S2P1 (diagram b) from the secondary interface.

variations is more complex with increasing the depth. There is a relatively small perturbation in the horizontal plane and a moderate velocity gradient in the vertical direction.
First, we give an example showing the primary reflected ray paths only (Fig. 15a). In the model parametrization the length of the cubic cell is 1.0 km and the number of secondary nodes (M) added onto each direction of the cell surface is 9. Two undulating interfaces are located at depths of 14 and 26 km, respectively. The sampled node spacing on the subsurface interfaces also has the same

spacing interval as the cell surface. The source was set at the top surface centre (X = 20.0, Y = 20.0 and Z = 0.0 km) of the model and 62 receivers were placed at the three sides of the top surface
boundaries with 2.0 km spacing interval. Two kinds of primary Pwave reflections (P1 P1 reflected phase from the upper subsurface interface and transmitted and reflected phase P1 P2 P2 P1 from the lower subsurface interface) were simulated.
Second, we show the ray paths for the same 3-D model, for
multiply transmitted, converted and reflections, but with horizontal

C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

210 C.-Y. Bai, X.-P. Tang and R. Zhao

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Figure 13. Ray paths of transmitted and converted primary reflected arrivals P1S2P3S3P2S1 (diagram a) and P1P2S3P3S2P1 (diagram b) from the third interface.
subsurface interface instead (Fig. 15b). The ray paths of reflected and converted phase P1 S1 from the upper horizontal subsurface interface and the multiple transmitted, converted and reflected phase P1 S2 P2 S2 P2 S1 from the lower horizontal subsurface interface are shown.

Figure 14. 3-D velocity model. Diagram (a) is a plan view of the velocity model in the model base (Z = –30.0 km). Dashed contours denote a velocity high, that is, a central high and decreasing values outward (indicated in figure); the continuous contours denote a velocity low, that is, a central low and increasing values outward. Diagram (b) shows the velocity model in X –Z cross-section (Y = 4.0 km).

5 ERROR ESTIMATIONS
The upper error bound with the one-stage MSPM approach in the computation of first arrivals has been estimated in 2-D/3-D velocity models (Bai et al. 2007). The main features of the multistage MSPM are the same as for the one-stage MSPM scheme. Therefore it can be expected that the multistage MSPM scheme would have the same order of error estimation in tracing multiple arrivals as is incurred in tracing first arrivals with the one-stage MSPM approach. Here, we give two examples of computing traveltime errors with the multistage MSPM scheme, determined by comparing numerical traveltimes against analytic solutions. First, we use the Model of Fig. 6 as a 2-D example, except here the layered velocity model is replaced by a homogenous velocity field and the three curved interfaces are replaced by three horizontal interfaces (reflection due to density contrast alone). The source was set at the top centre of the model and 51 receivers were located along the top surface at 2.0 km spacing. The relative errors (per cent) for the reflected traveltimes from the interfaces 1, 2 and 3 measured against the analytic solutions are depicted in Fig. 16. The second example is a 3-D one, based on the uniform Model with scale length of 100 km × 100 km × 40 km and two horizontal planes were located at the depth of 15 and 30 km, respectively. In model parametrization the cell size is 4 km × 4 km × 2 km and the number of secondary nodes added in each direction of the cell surface is 9. The source was set at one top centre of the model boundary (X = 50.0 km, Y = 0.0 km and Z = 0.0 km) and 441 receivers were uniformly located on the top surface of the model. The relative errors
C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Multiple arrival tracking with MSPM 211
Figure 15. (a) The 3-D ray paths of the transmissions and primary reflections. The dark lines are the ray paths of primary reflected phase P1P1 from the first undulated subsurface interface, and the thin lines are the ray paths of the multiple transmissions and primary reflections P1P2P2P1 from the second undulated subsurface interface. The source was set at the central top surface and 62 receivers were arranged along the three boundaries of the top surface with 2.0 km spacing interval. (b) The 3-D ray paths of the multiple transmissions, conversions and reflections. The dark lines are the ray paths of converted and reflected phase P1S1 from the first horizontal interface, and thin lines are the ray paths of the multiple transmitted, converted and reflected phase P1S2P2S2P2S1 from the second horizontal interface. The source and receiver layout are the same as in Fig. 15(a).
Figure 16. The relative traveltime error (per cent) versus distance for different depth reflectors measured against the exact (analytic) solutions (2-D case, see Model of Fig. 6). The dotted, thin and dark lines are the results for the P-wave reflected back from the first, second and third horizontal interfaces, respectively. C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

212 C.-Y. Bai, X.-P. Tang and R. Zhao

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

Figure 17. The relative traveltime error (per cent) versus distance for different depth reflectors measured against the exact (analytic) solutions in 3-D case. The left- and right-hand panels are the results for the P-wave reflected back from the first and second subsurface interfaces, respectively. The error magnitudes are indicated in the figure with different sizes.

(per cent) of the reflected traveltimes (reflections being due entirely to density change) from the horizontal subsurface interfaces 1 and 2, compared to the exact (analytic) reflected traveltimes, are depicted in Fig. 17. From Figs 16 and 17, we see that in both the 2-D and 3-D situations the traveltime errors are symmetric about the source and generally less than 0.125 per cent. In addition there is no obvious increase in error propagation with increasing depth to the reflector, which means that the numerical calculation is stable and robust. Note that the relative errors can be further reduced by decreasing the cell size or the secondary node spacing (see Fig. 2 and Table 1).
In practical applications, the problem is how to select a suitable cell size and number of secondary nodes. The influence of number of secondary nodes was investigated in the previous benchmark tests. Therefore, here we concentrate on the impacts of the cell size selection. In this test, using model of Fig. 6 as an example (2-D case), the source is placed in the upper left-hand corner and the 51 receivers are located along the top surface from positions X = 0.0 to 100.0 km, at 2 km spacing. To test the influence of the cell size, we keep the number of secondary nodes along each cell boundary constant (M = 9), and progressively halve the cell size, going from 10 to 0.625 km. In this design, we have five different lengths of cell (10.0, 5.0, 2.5, 1.25 and 0.625 km) to parametrize the velocity model. Taking the results of the smallest cell-size model as the reference value, we show in Fig. 18 the relative error at each receiver location for the three different curved depth reflectors (diagram a, for the first interface; diagram b, for the second interface and diagram c: for the third interface). From Fig. 18 we can draw the following conclusions: (1) there is not much difference in the error estimations for the different depth reflectors. In other words, the error propagation is well constrained with distance; (2) apart from the largest cell size (10.0 km), the relative errors against the referenced values of fine cell-size model can be controlled within 0.1 per cent, independent of interface depth; (3) there exists a certain range for cell size selection. Reducing the cell size further does not sufficiently improve the accuracy as one might expect.
Such conclusions coincide with our previous results for one-stage MSPM in the calculation of the first arrivals (Bai et al. 2007). As

Figure 18. The variations of relative traveltime error versus distance with different cell-size models against the reference values (diagrams a, b and c: results of the P-wave reflected from the first, second and third interfaces, respectively).
shown previously, and in this study, the preferred approach is to choose a suitable cell size according to a specified velocity model and then add sufficient secondary nodes along each cell boundary (2-D case) or along each cell surface (3-D case) to guarantee a satisfactory accuracy in the computation. This applies to both the one-stage and the multistage MSPM.

C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

6 DISCUSSIONS AND CONCLUSIONS
We have presented a multistage MSPM scheme for tracing any kind of combination of multiply transmitted, reflected and mode converted arrivals in complex 2-D/3-D layered media and conducted a comparison with the multistage FMM approach in 2-D case. The desirable features when the multistage scheme is incorporated with the MSPM approach include (1) stability and a global solution as in the one-stage MSPM approach, (2) the ability to track multiple arrivals composed of any number of transmitted, reflected (or refracted) and converted branches, (3) a higher accuracy can be expected compared with the FMM algorithm for tracking the multiple arrivals and (4) less CPU time is consumed than with the multistage FMM approach. Application of the multistage MSPM scheme to highly heterogeneous media (such as the Marmousi model) having severe variations in wavespeed and interface curvature supports these conclusions.
Numerical results suggest that a suitable cell size and a sufficient number of secondary node supplements can retain a satisfactory accuracy. Although later arrivals (such as swallowtails) caused by continuous variations in wavespeed cannot be traced with the current multistage MSPM scheme, phases generated by velocity interfaces form a major class of multiple arrivals. As such, the new approach may be applied in a wide range of areas including multiple traveltime tomography (or multiple reflection tomography), the processing of coincident reflection events in seismic exploration data, and precise earthquake location. Future work should include how to incorporate reflection/transmission coefficients and attenuation for seismic amplitude studies.
ACKNOWLEDGMENTS
This research work was supported jointly by the National Science Foundation of China (NSFC, project no. 40534021 and 40774020) and the Key Research Program of the Ministry of Education of China (MEC, project no. 107137). The authors would like to thank Dr Nicholas Rawlinson in RSES of Australian National University for generously offering their 2-D codes of multistage FMM approach and detailed readme text for how to use them in the PC environment and also two anomalous reviewers whose suggestions improved the revised manuscript substantially. Additionally, the authors would like to thank Dr Laura Codben in SPICE, who was helpful in polishing the English presentation of the text.
REFERENCES
Afnimar & Koketsu, K., 2000. Finite difference traveltime calculation for head waves travelling along an irregular interface, Geophys. J. Int., 143, 729–734.
Bai, C.-Y. & Greenhalgh, S.A., 2005a. 3-D non-linear travel time tomography: imaging high contrast velocity anomalies. Pure appl. Geophys., 162, 2029–2049.
Bai, C.-Y. & Greenhalgh, S.A., 2005b. 3-D multi-step travel time tomography: imaging the local, deep velocity structure of Rabaul volcano, Papua New Guinea, Phys. Earth planet. Int., 151, 259–275.
Bai, C.-Y., Greenhalgh, S & Zhou, B., 2007. 3-D ray tracing with a modified shortest path method, Geophysics, 72, T27–T36.
Bai, C.-Y., Peng, J.-B. & Grenhalgh, S., 2008. 3-D P-wave velocity structure of the crust beneath the Loess Plateau and surrounding regions of China, Tectonophysics, 460, 278–287.
Cao, S. & Greenhalgh, S., 1993. Calculation of the seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm, Geophys. J. Int., 114, 593–600.
C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS

Multiple arrival tracking with MSPM 213
Cao, S. & Greenhalgh, S., 1994. Finite difference solution of the eikonal equation using an efficient, first arrival wavefront tracing scheme, Geophysics, 59, 632–643.
Ceˇrveny , V., 1987. Ray tracing algorithms in three-dimensional laterally varying layered structures, in Seismic Tomography: With Applications in Global Seismology and Exploration Geophysics, pp. 99–133, ed. Nolet, G., D. Reidel Publishing Company, Amsterdam, The Netherlands.
Ceˇrveny’, V., 2001, Seismic Ray Theory, 1st edn, Cambridge University Press, Cambridge, UK.
Cheng, N. & House, L., 1996. Minimum traveltime calculations in 3-D graph theory, Geophysics, 61, 1895–1898.
Coman, R. & Gajewski, D., 2005. Traveltime computation by wavefrontorientated ray tracing, Geophys. Prospect., 53, 23–26.
De Kool, M., Rawlinson, N. & Sambridge, M., 2006. A practical gridbased method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media, Geophys. J. Int., 167, 253–270.
Di Stefano R., Aldersons, F., Kissling, E., Baccheschi, P., Chiarabba, C. & Giardini, D., 2006. Automatic seismic phase picking and consistent observation error assessment: application to the Italian seismicity, Geophys. J. Int., 165, 121–134.
Fischer, R. & Lees, J.M., 1993. Shortest path ray tracing with sparse graphs, Geophysics, 58, 987–996.
Fomel, S. & Sethian, J.A., 2002. Fast-phase space computation of multiple arrivals, Proc. Natl. Acad. Sci. U.S.A., 99, 7329–7334.
Gruber, T. & Greenhalgh, S., 1998. Short note: precision analysis of firstbreak times in grid models, Geophysics, 63, 1062–1065.
Harten, A., Engquist, B., Osher, S. & Chakravarthy, S.R., 1987. Uniformly high order accurate essentially non-oscillatory schemes, III, J. Computat. Phys., 71, 231–303.
Hauser, J., Sambridge, M. & Rawlinson, N., 2008. Multiarrival wavefront tracking and its applications, Geochem. Geophys. Geosyst., 9, Q11001, doi:10.1029/2008Gc002069.
Hole, J.A. & Zelt, B.C., 1995. 3-D finite-difference reflection travel times, Geophys. J. Int., 121, 427–434.
Jiang, G.S. & Peng, D.P., 2000. Weighted ENO schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput., 21, 2126–2143.
Jiang, G.S. & Shu, C.W., 1996. Efficient implementation of weighted ENO schemes, J. Computat. Phys., 126, 202–228.
Julian, B.R. & Gubbins, D., 1977. Three-dimensional seismic ray tracing, J. Geophys., 43, 95–113.
Kennett, B.L.N. & Engdahl, E.R., 1991. Traveltimes for global earthquake location and phase associations, Geophys. J. Int., 105, 429–465.
Kim, S., 2002. 3-D eikonal solvers: first-arrival traveltimes, Geophysics, 67, 1225–1231.
Klimesˇ, L. & Kvasnicˇka, M., 1994. 3-D network ray tracing, Geophys. J. Int., 116, 726–738.
Lambare´, G., Virieux, J., Madariaga, R. & Jin, S., 1992. Iterative asymptotic inversion in the accoustic approximation, Geophysics, 57, 1138– 1154.
Liu, X., Osher, S. & Chan, T., 1994. Weighted essentially non-oscillatory schemes, J. Computat. Phys., 115, 200–212.
Moser, T.J., 1991. Shortest path calculation of seismic rays, Geophysics, 56, 59–67.
Nakanishi, I. & Yamaguchi, K., 1986. A numerical experiment on nonlinear image reconstruction from first-arrival times for twodimensional island arc structure, J. Phys. Earth, 34, 195–201.
Osher, S., Cheng, L.T., Kang, M., Shim, Y. & Tsai, Y.H., 2002. Geometric optics in a phase-space-based level set and Eulerian framework, J. Computat. Phys., 179, 622–648.
Podvin, P. & Lecomte, I., 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools, Geophys. J. Int., 105, 271–284.
Qian, J. & Symes, W.W., 2002. An adaptive finite-difference method for traveltimes and amplitudes, Geophysics, 67, 167–176.
Rawlinson, N. & Sambridge, M., 2004. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method, Geophysics, 69, 1338–1350.

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

214 C.-Y. Bai, X.-P. Tang and R. Zhao
Reading A.M., Mao W.-J. & Gubbins D., 2008. Polarization filtering for automatic picking of seismic data and improved converted phase detection. Geophys. J. Int., 147, 227–234.
Riahi, M.A. & Juhlin, C., 1994. 3-D interpretation of reflected arrival times by finite-difference techniques, Geophysics, 59, 844–849.
Schneider, Jr. W.A., 1995. Robust and efficient upwind finite-difference traveltime calculations in three dimensions, Geophysics, 60, 1108– 1117.
Sethian, J.A., 1996. A fast marching level set method for monotonically advancing fronts, Proc. Nat. Acad. Sci., 93, 1591–1595.
Shu, C. & Osher, S., 1989. Efficent implementation of essentially nonoscillatory shock-capturing schemes, II, J. Computat. Phys., 83, 32– 78.
Um, J. & Thurber, C., 1987. A fast algorithm for two-point seismic ray tracing, Bull. seism. Soc. Am., 77, 972–986.

Van Avendonk, H.J.A., Harding, A.J., Orcutt, J.A. & Holbrook, W.S., 2001. Hybrid shortest path and ray bending method for traveltime and ray path calculation, Geophysics, 60, 648–653.
van Trier, J. & Symes, W.W., 1991. Upwind finite-difference calculation of traveltimes, Geophysics, 56, 812–821.
Vidale, J.E., 1988. Finite-difference calculations of traveltimes, Bull. seism. Soc. Am., 78, 2062–2076.
Vidale, J.E., 1990. Finite-difference calculations of traveltimes in three dimensions, Geophysics, 55, 521–526.
Vinje, V., Iversen, E. & Gjoystdal, H., 1993. Travel-time and amplitude estimation using wave-front construction, Geophysics, 58, 1157–1166.
Zhang, J. & Toksoz, M.N., 1998. Nonlinear refraction traveltime tomography, Geophysics, 63, 1726–1737.
Zhao, D., Tode, S. & Lei, J., 2005. Local earthquake reflection tomography of the Lander aftershock area, Earth planet. Sci. Lett., 235, 623–631.

Downloaded from https://academic.oup.com/gji/article/179/1/201/735460 by guest on 02 March 2024

C 2009 The Authors, GJI, 179, 201–214 Journal compilation C 2009 RAS