Coupling spectral elements and modes in a spherical Earth: an extension to the ‘sandwich’ case

SUMMARY We present an extension to the coupling scheme of the spectral element method (SEM) with a normal-mode solution in spherical geometry. This extension allows us to consider a thin spherical shell of spectral elements between two modal solutions above and below. The SEM is based on a high-order variational formulation in space and a second-order explicit scheme in time. It combines the geometrical ﬂexibility of the classical ﬁnite-element method with the exponential convergence rate associated with spectral techniques. In the inner sphere and outer (cid:2)(cid:2) -layer model based on the tomographic model SAW24b16 is presented up to a corner frequency of 1/12 s. The comparison with data shows surprisinglygoodresultsforthe3-Dmodelevenwhentheobservedwaveformamplitudesdiffersigniﬁcantlyfromthosepredictedinthesphericallysymmetricreferencemodel(PREM).

Peak-to-peak lateral variations of up to 10 per cent in S velocity occurring over several hundred km have been found, both on the border of the African Plume (e.g. Ritsema et al. 1998;Wen & Helmberger 1998a;Ni & Helmberger 2001) and of the Pacific Plume (Bréger & Romanowicz 1998;Bréger et al. 1998). Such strong variations cannot be interpreted in terms of thermal variations alone. Strong localized variations in P velocity have also been inferred from the study of precursors to PKP (Vidale & Hedlin 1998;Wen & Helmberger 1998a) and in S -velocity anisotropy from the study of SV diff (Vinnik et al. 1998).
Yet, present global waveform modelling approaches rely heavily on assumptions of weak heterogeneity. While the forward modelling of traveltimes of body wave phases sensitive to the base of the mantle, using standard ray methods, provides helpful insights regarding the character of lateral heterogeneity, much information is yet to be gained from the analysis of waveforms. For this purpose, adequate modelling tools need to be applied.
Given the strong heterogeneity found in the two boundary layers of the mantle, appropriate tools are needed that will handle waveform modelling of: (1) the propagation of seismic body and surface waves in 3-D models with strong lateral variations and in spherical geometry and (2) the diffracted waves along the coremantle boundary. Diffracted waves cannot be handled by ray-based methods. On the other hand, perturbation methods based on a normal-mode formalism are well adapted to the spherical geometry, can handle diffracted waves and allow the computation of Fréchet kernels for inversion (e.g. Lognonné & Romanowicz 1990;Li & Romanowicz 1995Clévédé & Logonnné 1996;Dahlen et al. 2000). However, the strength of the target lateral heterogeneity, which would require pushing the perturbation development to rather large orders beyond the Born approximation, and the relatively short period of the waves considered (30 s or less), makes this approach as yet rather impractical for waveform modelling, especially in the P -SV case (many modes to couple). Another approach that has been proposed for whole Earth heterogeneous models is based on the direct solution method (DSM, Geller & Ohminato 1994;Geller & Takeuchi 1995). This method, based on the weak form of the equations of motions, allows one to compute partial derivatives of seismograms, and is therefore well adapted for inversion. Unfortunately, this method is currently only available for axisymmetric models (e.g. Cummins et al. 1997) and uses the Born approximation for models without symmetry (Takeuchi & Geller 2000).
There have been some successes in modelling D sensitive phases using hybrid codes, in which the computation in D is performed using a Cartesian finite-difference (FD) scheme, while outside of the deep mantle, standard 1-D ray methods such as WKBJ (Chapman 1978), Kirchoff (Stead & Helmberger 1988) or generalized rays (Helmberger 1983), are applied to perform the ray tracing. The FD part is computationally intensive and the hybrid approach reduces the computation time outside of the target heterogeneous region. Wen & Helmberger (1998b,a) successfully implemented such a dual scheme, and in particular modelled the effect of the Ultra Low Velocity Zone (ULVZ) on the PKP and SKS + SPdKS waveforms. Because the disturbance due to the ULVZ is limited to a small region of D , this type of approach is, from a computational perspective, particularly efficient. It is unfortunately not appropriate for modelling S diff waves. First, S diff can diffract over more than 20 • epicentral distance, and it is not possible to adequately simulate diffraction over a curved CMB in the Cartesian FD box.
Over the last few years, much progress has been made in the development of numerical methods adapted to spherical geometry and able to compute waves emanating from a realistic seismic source, reaching, within reasonable computational time, periods of interest for teleseismic studies, making no assumptions on the strength of velocity contrasts, and able to handle interface waves and interface topography.
Among the possible numerical methods able to solve the wave equation in general Earth models, the spectral element method (SEM) has been shown to be particulary efficient and accurate. The SEM has been introduced in computational fluid dynamics (Patera 1984;Maday & Patera 1989) and applied more recently to the 3-D elastic equation (Faccioli et al. 1997;Komatitsch & Vilotte 1998). This method combines the geometrical flexibility of conventional finite-element methods with the exponential convergence rate associated with spectral techniques, and suffers from minimal numerical dispersion and diffusion. The extension to spherical geometry has been introduced by Chaljub (2000) and Chaljub et al. (2003), developing a mesh of the sphere with deformed cubes named the 'cubic sphere' starting from the work of Sadourny (1972) and Ronchi et al. (1996), and allowing non-conforming interfaces using the mortar method (Bernardi et al. 1994). The effects of anisotropy, attenuation (Komatitsch & Tromp 1999), rotation and gravity (Komatitsch & Tromp 2002) have also been introduced. In spite of all the qualities of the method, the main drawback is still the numerical cost. The method can address corner frequencies between 1/20 and 1/15 Hz but only with huge amounts of memory (of the order of 100 Gb) and with a large CPU time, making it still impractical to test numerous models, as one would want to do in a forward modelling approach, or compute the wavefield for hundreds of sources in an iterative global inversion scheme.
A solution to that problem, allowing higher-frequency simulations using smaller machines (smaller memory and less CPU time), has been introduced with the coupling of the SEM with the modal solution (Capdeville 2000;Capdeville et al. 2002Capdeville et al. , 2003. The idea of this method is to limit the use of the expensive SEM in regions of the Earth that, depending on the problem studied, include 3-D heterogeneity, and to use the cheaper modal solution in regions that can be assumed to be spherically symmetric. The coupling between the SEM, expressed in the space and time domain, and the modal solution, expressed in the frequency and wavenumber domain, is not trivial and requires some original solutions that are explained in detail in Capdeville et al. (2003)(hereafter referred to as 'Paper 1'). In the spectral element method, the coupling is introduced via a dynamic coupling operator, a Dirichlet-to-Neumann (DtN) boundary operator. The operator can be explicitly constructed in the frequency domain and in the generalized spherical harmonics basis, using classical modal solution techniques.
In Paper 1, only one coupling interface was allowed with an external shell of spectral elements over an inner sphere in which the modal solution is found. The typical application of such a partitioning is for a heterogeneous crust over a spherically symmetric mantle and core. Here we present an extension of the coupled method to the case of a thin spherical shell of spectral elements 'sandwiched' between two modal solutions. A target application of this new development is the study of the D layer, which is strongly heterogeneous, as discussed above. In that case, a layer with 3-D structure at the bottom of the mantle can be used between two spherically symmetric models. This particular coupling type is not trivial, due to the fact that the source and the receivers are in the modal solution domain, as we will show here. We finally present and discuss an example of simulations, with a 3-D model in the D layer obtained using the S tomographic model SAW24b16 (Mégnin & Romanowicz 2000) C 2003 RAS, GJI, 154, 44-57 Figure 1. The Earth domain Ω (left) is divided into three parts (right), an external shell Ω M2 , an internal shell Ω S and an internal sphere Ω M1 separated by two spherical boundaries Γ 2 and Γ 1 . In this sketch, Γ 1 is located on the core-mantle boundary and Γ 2 in the lower mantle. We assume that lateral heterogeneities of the Earth model are only present in Ω S , so that the SEM needs only to be used in that domain and modal solutions in Ω M1 and Ω M2 . and PREM (Dziewonski & Anderson 1981) and some comparisons with data.

P RO B L E M S TAT E M E N T A N D M E T H O D P R I N C I P L E
We consider a non-rotating Earth Ω of radius r Ω . The equation of motion to be solved in Ω is where ρ is the density, u is the displacement field,ü is its second partial derivative with respect to time t, H is the elasto-gravity operator for a non-rotating spherical Earth (e.g. Valette 1986;Woodhouse & Dahlen 1978) and f the generalized body force due to the earthquake. We assume a free surface boundary condition ∂Ω, and an initial state of the form u(r, 0) =u(r, 0) = 0. We divide the Earth into three parts: an external shell Ω M2 , an internal sphere Ω M1 and an internal shell Ω S sandwiched between them. We name Γ 1 the spherical boundary between the domains Ω M1 and Ω S , and Γ 2 the spherical boundary between the domains Ω S and Ω M2 (see Fig. 1). Note that if Ω M2 is reduced to 0, we are in the same situation as presented in Paper 1, that is one external shell over an inner sphere. We assume that all lateral heterogeneities of the Earth model are localized in the inner shell Ω S , and that the Earth model within Ω M1 and Ω M2 , as well as on Γ 1 and Γ 2 , is spherically symmetric. Depending on the Earth model considered, the radius r Γ i of interface Γ i can be set anywhere between the core-mantle boundary radius and r Ω . The basic idea of the method is to use the spectral element method in the heterogeneous part, that is Ω S , and a modal solution in Ω M1 and Ω M2 . The latter is well known when the model properties are only varying radially. The point is that, on the one hand, the spectral element method is very well adapted to a general 3-D medium but is time and memory consuming from a numerical point of view. On the other hand, the modal solution has a very low numerical cost in spherically symmetric media. Combining these two methods, we expect to optimize the numerical cost of wave propagation in Earth models where, for example, we wish to focus the investigation on lateral heterogeneities in a given depth range of the mantle, such as D and its vicinity.

Variational formulation
We first solve the wave equation using the SEM in Ω S . The coupling with the modal solution in Ω M1 and Ω M2 will then come naturally. The SEM is a finite-element method that solves the equation of motion in its variational form, which is the integral form of eq. (1).
In this problem formulation, we seek a solution in V, the set of square-integrable functions with square-integrable generalized first derivatives over Ω S . The problem to be solved is: find u(·, t) ∈ V, such that ∀t ∈ I = [0, T ], the time duration of the simulation, and ∀w ∈ V is the classical L 2 inner product, the symmetric elasto-gravity bilinear form a (·,·) expression can be found in Paper 1 and where T Γ i is the traction on the spherical interface Γ i . To solve eq.
(2), we need to know T Γ i . As will be shown in Section 2.3, T Γ i can be computed as a function of the incident displacement u Γ i on each coupling interface and this is where the coupling is performed.

The spectral element approximation
The basic principle of the SEM, which is very close to the classical finite-element method, is to solve eq. (2) using a high-degree polynomial approximation by elements of functions in V space. In this method, elements have to be deformed cubes, so a cubic meshing of the spherical shell has to be found. The 'cubic sphere' proposed by Sadourny (1972) and further extended by Ronchi et al. (1996) allows such a meshing of a spherical surface by decomposing it into six regions of identical shape, which can be mapped on to a cube face. To obtain the meshing of a spherical shell, spherical surfaces are connected radially (see Fig. 2), where non-conforming interfaces are allowed . Each numerical integration of eq. (2) is performed using the Gauss-Lobatto-Legendre (GLL) quadrature in each Cartesian direction. The polynomial basis is built using the Lagrange polynomial associated with GLL points. A detailed description of the spectral element method applied to the wave equation can be found in Komatitsch & Vilotte (1998) and Chaljub et al. (2003). In this paper, the anelasticity of the medium is taken into account in the SEM following the scheme presented in Komatitsch & Tromp (1999).

The Dirichlet to Neumann operator
In this part, we recall results obtained in Paper 1 ignoring for the moment that the source and the receivers will be in the upper modal solution domain, Ω M2 . This aspect will be examined in the next section. The continuity of traction and the continuity of displacement, or the normal displacement (depending on whether Γ is a solid-solid or a solid-fluid interface), through each interface Γ i , i = 1, 2, have to be assured. Assuming that the solution to the wave equation (1), without the right-hand side f, is known in Ω Mi , if a boundary condition in displacement is imposed on the interface Γ i , the stress field can be computed everywhere in Ω Mi and, in particular, the traction on Γ i is known. Using the continuity relations of displacement and traction through Γ i , we are then able to construct an operator A i that, for a given displacement u Γ i on Γ i , returns the corresponding traction that Ω Mi applies on Ω S : for a solid-solid interface and for a solid-fluid interface, where n is the normal outward unit vector to the surface Γ i . Form the point of view of Ω S , the boundary condition on Γ i is a Neumann condition (condition in traction) that depends on a Dirichlet condition (condition in displacement), therefore the operator A i is named a Dirichlet-to-Neumann operator. This operator allows us to compute eq. (3) knowing u Γ i , as for a classical absorbing boundary problem (e.g. Givoli & Keller 1990;Grote & Keller 1995;Sánchez-Sesma & Vai 1998).
As shown in Paper 1, the DtN operator is built in the frequencygeneralized spherical harmonic domain (Phinney & Burridge 1973), in which the solutions of eq. (1) without the right-hand side term are well known in spherically symmetric models. For each angular order , an operator in frequency A (ω) is found which, due to the spherical symmetry of the problem, does not depend on the azimuthal order (m). This operator is not defined for a discrete set of frequencies d , which correspond to the eigenfrequencies of Ω Mi for the homogeneous Dirichlet problem (no displacement at the boundary Γ i ). Note that if one needs to include attenuation in the spherically symmetric part of the model, anelasticity is simply approximated by introducing a complex part in each eigenfrequency  of d , as is classically done in normal-mode problems (e.g. Takeuchi & Saito 1972). Because the SEM is a time-space method, the DtN operator has to be computed in this domain. The first and most difficult step is to compute A in the time domain for each . First, because of its singularities at each frequency of d , this operation cannot be performed by a traditional fast Fourier transform (FFT). To circumvent this problem, the continuous spectrum, on which a classical FFT can be performed, is separated form the discrete spectrum. The Fourier transform of the discrete part of the operator can be obtained using the Cauchy theorem, which finally allows us to obtain the DtN operator in the time domain. The DtN operator is naturally causal, but in order to be compatible with the SEM time evolution scheme, it has to be numerically causal. By numerically causal we mean that the DtN operator in time that is obtained from the DtN operator in frequency, using the numerical process previously described, must be equal to zero before t = 0. If this were not the case, u Γ i (t) would be required at time steps in advance of the current time step, in order to compute T Γ i (t) correctly at the current time step, which is obviously not possible. Unfortunately, the frequency window used for the Fourier transform is not infinite (the maximum frequency possible is the Nyquist frequency), and such a filtering is not causal. In order to circumvent this second problem, a regularized DtN operator, A r , is used, where A r (ω) is close to zero for the high frequencies. On such a regularized DtN operator, the frequency window filtering has no effect and the causality is preserved numerically. A r is obtained by subtracting from A an asymptotic operator C valid for the high frequency of the DtN operator, which can be obtained analytically (see Paper 1). With such a regularization, the expression of the traction in the generalized spherical harmonic domain is where * is the time convolution. As A r is causal, A r * u Γ i ,m can be computed numerically without any problem and it can be shown that C * u Γ i , ,m can be computed analytically. Mathematically speaking, A is not a bounded operator and therefore, the numerical Fourier transform cannot be performed on it properly. The regularized operator A r is bounded (it goes down to zero for high frequency), and therefore the numerical Fourier transform can be performed on it properly, which solves the problem. After having obtained T Γ i , ,m in time, the second step is to obtain it in space. To do so, a backward Legendre transform, that is the summation over and m of coefficients T Γ i , ,m on the generalized spherical harmonic basis, has to be performed. The summation over has to be numerically truncated after an max that does not affect the coupling process (the summation over m is naturally truncated as, for a given , m must lie between − and ). To evaluate this corner angular order max , the dispersion curve of the surface waves of the inner sphere for the homogeneous Dirichlet boundary condition problem can be used. As a matter of fact, with such a curve, for a given maximum frequency of the source, a maximum angular order can be found. This maximum angular order corresponds to the maximum angular order that a wave would have in the inner sphere in the far field of the source. Multiplying this angular order with a 'safety' coefficient (Γ i can be in the near field of the source and the medium close to the interface in Ω S can be strongly heterogeneous), let us say 2, a very good coupling is obtained. Finally, the traction expression which is to be used to evaluate eq. (3) is where Y ,m is the generalized spherical harmonic basis.

Particularity of the coupling due to the 'sandwich'
The introduction of the 'sandwich' geometry has two practical consequences that have to be treated: the source and the receivers are in the upper modal solution domain. The resolution and the construction of the DtN operator remains unchanged in the inner sphere Ω M1 , but in the outer shell Ω M2 , we must take into account the right-hand sidefin the resolution of eq. (1). This is performed by adding a particular solution of eq. (1) to the general one which gives the following DtN relation: in the solid-solid case (the solid-liquid case is similar), where B is the particular solution term due to the presence of the source in Ω M2 . B is computed using the modal solution of Ω M2 where the boundary condition on Γ 2 is chosen as a homogeneous Dirichlet condition (no displacement) for practical reasons (see the Appendix for details). As receivers are on the free surface of the Earth, they are located in the modal solution domain Ω M2 . We must therefore use the modal solution to obtain the displacement at the surface. This is performed using an operator P similar to the DtN operator that continues the SEM domain solution on the coupling interface Γ 2 up to the free surface in the modal solution domain: where u M2 is the displacement in Ω M2 , B d is a term similar to B but in displacement and at the free surface. All operators A 2 , B, B d and P have a discrete spectrum that is a subset of the spectrum of the spherical shell Ω M2 with a free surface condition at the surface and a homogeneous Dirichlet condition on Γ 2 but they are not all exactly the same. Indeed, some eigenfrequencies corresponding to surface waves do not contribute to the DtN after a certain frequency depending on the depth of Γ 2 , and are therefore not present in the spectrum of A 2 , B and P. However, they are present in the spectrum of B d . Construction of operators B and P is detailed in the Appendix.

VA L I D AT I O N T E S T S
Before presenting a test of the 'sandwich' coupling, we first present a validation test for a diffracted wave on the coupling interface. A diffracted wave exactly on the coupling interface is a very difficult case for this method because such a wave stays on the DtN boundary for a long time and is therefore very sensitive to any error that occurs during the coupling process (an even worse case, i.e. Stoneley waves on the coupling boundary, is discussed in Capdeville 2000). Furthermore, this kind of wave (S diff and P diff ) is widely used to study the D layer, and because we do not consider any heterogeneity in the outer core, the coupling interface Γ 2 will be set at the CMB and therefore diffracted waves will propagate on the DtN boundary. It is therefore crucial to test the accuracy of simulations in this particular case.
To do so, we use a simple model that is a homogeneous spherical shell (external radius 6371 km) over a liquid inner sphere (radius 2871 km). The external shell S -wave velocity (β) is 6 km s −1 , the P-wave velocity (α) is 8 km s −1 , the inner sphere P wave is set to 4 km s −1 to create a shadow area for P waves so that P diff and PKP waves do not mix (Fig. 3). The density is everywhere 3000 kg m −3 . The source, located at a depth of 1048 km, is an explosion C 2003 RAS, GJI, 154, 44-57 Downloaded from https://academic.oup.com/gji/article-abstract/154/1/44/604492 by guest on 29 July 2018 (this implies that no SH diff waves will be generated, but SH diff is not a difficult case because for this wave, the boundary is only a free surface), and the corner frequency is 1/125 Hz. Fig. 4 shows the spectral element mesh used for that test. We compare synthetics obtained with the coupled method with normal-mode synthetics at four different epicentral distances on Fig. 5. The residuals, amplified 10 times, show a very good agreement, which validates the coupling in that case and demonstrates that diffracted waves propagating along the coupling interface are computed with a satisfactory accuracy. More tests of the DtN but non-specific to the sandwich coupling can be found in Capdeville (2000) or Capdeville et al. (2003).
To validate the sandwich coupling, we perform a test in a fully homogeneous sphere with a very deep source. The test is once again unrealistic geophysically speaking, but has every difficulty, and even more (the source is not usually so deep), of a realistic case for the coupling. Furthermore, the normal-mode solution is, in that case, quasi-analytic and therefore suitable. The source is very deep in order to minimize surface waves that could hide problems at the coupling interface. The elastic properties of the sphere are the same as those of the external shell in the previous test. The different radii are r Ω = 6371 km, r Γ 2 = 3810.5 km and r Γ 1 = 2560.5 km. The spectral element mesh is exactly the same as that in Fig. 4, but, because of geometrical effects (radii are smaller), the maximum corner frequency can be higher (1/55 Hz). The source is a strike-slip earthquake at 1272 km depth. In Fig. 6, we present the contribution of the two terms of eq. (9) to the actual seismogram. The term B d represents the source contribution in a spherical shell with a rigid boundary condition at the bottom. The reflection of the inner interface can be clearly seen, especially on the transverse component. The term P * u M2 Γ 2 is the contribution of the DtN operator to the final seismogram. This contribution cancels all the reflections at the inner interface of the term B d , to finally obtain a very good match with the normal-mode summation reference solution. Once again the accuracy is satisfactory.
Finally, we show in Fig. 7 a comparison with the normal-mode solution in PREM. The configuration is similar to that in the previous test, but the source is now 600 km deep and r Γ 1 has been set to 3480 km to match the core-mantle boundary. Once again, residuals for both vertical and transverse components show a very good agreement.

A N E X A M P L E O F A P P L I C AT I O N : T H E D L AY E R
We perform a simulation of one deep event (621 km deep) of magnitude 6.8 (1997 September 4, Fiji) using PREM in the top shell, and the 3-D degree-24 SH model SAW24b16 (Mégnin & Romanowicz 2000), in the 370 km above the CMB. To obtain the P -wave velocity and density heterogeneities, we use simple linear relations, δρ = 0.4δβ and δα = 0.25δβ. The source mechanism is obtained from the Harvard CMT catalogue and we have chosen The spectral element mesh used has about 12 000 elements and 4 × 10 6 integration points (see Fig. 9), which allows a minimum corner period of 12 s.
We first present synthetics in PREM and SAW24b16 at two stations (INK and LMN) in Fig. 10. The effect on P waves (vertical component) is weak at this frequency, but the effect on S diff (transverse component) is noticeable, especially for the amplitude at LMN. This effect should be observable on the corresponding data.
Our next step is to compare data with synthetics produced in the 3-D model and in PREM for reference. To do so, the station response and a crude ellipticity correction have been applied to synthetics. Fig. 11 presents such a comparison for the 10 stations considered. These 10 stations have been chosen for their representativity when comparing data and synthetics. Two phases are shown, ScS or S diff and sScS or sS diff (depending on the epicentral distance). The first observation is that, for most of the stations and to first order, the 3-D model does a better job than PREM, both on the time delay and on the amplitude. This is especially true for stations such as LMN and BOSA at large epicentral distance, where the effect on the amplitude is the strongest. Note that slow regions are systematically associated with higher amplitude than PREM (e.g. YKW3, LMN) and fast regions with smaller amplitude than PREM (BOSA). The fact that a tomographic model gives a good result on the amplitude was not obvious a priori and is a good surprise. However, the 3-D models explain only the first-order features of the data, and not at all stations. The time-shift due to the 3-D model is sometimes too strong (e.g. WMQ) and sometimes both the amplitude and the phase are poorly explained (e.g. BW06 or TLY). It shows that interesting work still remains to better explain observed diffracted waveforms, even at the relatively long periods considered. The second phase, sS diff , is poorly modelled. This can be explained by the fact that this wave spends a significant amount of time in the strongly heterogeneous upper mantle near the source, and this heterogeneity is not accounted for in the model. This shows one of the limitations of our approach linked to the fact that we do not take into account heterogeneities anywhere other than in and above the D layer. Nothing prevents us, however, from progressively incorporating heterogeneity at different levels in the mantle, and in particular, from considering, in the future, two or more shells of strong 3-D structure, as needed.
Finally, in Fig. 12 we plot time arrivals of ScS or S diff phases computed by linearized ray tracing and by the coupled method using waveform cross-correlation for a large number of stations. It shows that in most cases ray tracing and the coupled method have time residuals of the same sign, but with significant differences of absolute value. The general trend is that linearized ray theory overestimates the time residual, which is coherent with the wave front healing phenomenon (Hung et al. 2001). A more extensive discussion of this type of comparison will be given in a forthcoming paper. This particular simulation was performed on 64 processors of the IBM machine of the NERSC, it required about 13 Gb of memory and lasted approximatively 20 h.

D I S C U S S I O N A N D C O N C L U S I O N S
We have presented an extension to the coupled spectral elements/modes method, which allows us to consider a thin spherical shell of spectral elements 'sandwiched' between two modal solutions. This extension provides a way to obtain relatively highfrequency seismograms at reasonable computation cost to study 3-D structure in specific shells of the Earth. The accuracy of the method is checked against normal-mode summation in simple models and shows satisfactory precision. An important application, as shown here, is the study of the D layer and its vicinity, where we can reach corner frequencies under 10 s on moderately large parallel computers (typically 64 processors and under 20 Gb of memory). Using this tool, we hope to provide strong constrains on the 3-D heterogeneity in and above D , in well-sampled regions of geodynamical interest, such as in the region of the Pacific superplume.
The comparison of observed S diffracted seismograms for paths sampling D across the Pacific, with synthetics computed in an existing tomographic model in which heterogeneity has been restricted to the bottom 370 km of the mantle shows surprisingly good agreement, not only in phase, but also in amplitude (in contrast with PREM synthetics), at least down to a 12 s corner frequency. This indicates that 3-D effects not accounted for by the theoretical approximations used in the construction of model SAW24B16 are not systematically dominant. Notable differences remain, and will be investigated further. The main limitation of the approach presented in this paper is of course the fact that the model is not 3-D everywhere, but this is the price to pay to be able to reach interesting frequencies and not be too restricted in the number of trial models to run.
Another interesting target of such a method is the inner core. In that case, the SEM would be used only in the inner core and the modal solution everywhere else. Corner frequencies of 5 s should be within range with the same type of machine.
Among the possible future developments, it is possible to introduce some 3-D structure in the modal part, such as ellipticity, using modal perturbations. However, it will not be possible to include general 3-D models without falling again in the classical difficulties of normal-mode perturbation techniques.

A C K N O W L E D G M E N T S
We thank Dimitri Komatitsch and Jeroen Tromp for kindly providing us with their attenuation scheme and code for the SEM. This work was partly supported by the Miller institute of Berkeley. and by NSF grant EAR-0106000. It is BSL contribution no 0303. All the simulations have been performed on 'seaborg', the IBM RS6000 supercomputer of the NERSC.