Summary
We present an azimuthally anisotropic 3-D shear-wave speed model of the Australian upper mantle obtained from the dispersion of fundamental and higher modes of Rayleigh waves. We compare two tomographic techniques to map path-average earth models into a 3-D model for heterogeneity and azimuthal anisotropy. Method I uses a rectangular surface cell parametrization and depth basis functions that represent independently constrained estimates of radial earth structure. It performs an iterative inversion with norm damping and gradient regularization. Method II uses a direct inversion of individual depth layers constrained by Bayesian assumptions about the model covariance. We recall that Bayesian inversions and discrete regularization approaches are theoretically equivalent, and with a synthetic example we show that they can give similar results. The model we present here uses the discrete regularized inversion of independent path constraints of Method I, on an equal-area grid. With the exception of westernmost Australia, we can retrieve structure on length scales of about 250 km laterally and 50 km in the radial direction, to within 0.8 per cent for the velocity, 20 per cent for the anisotropic magnitude and 20° for its direction. On length scales of 1000 km and longer, down to about 200 km, there is a good correlation between velocity heterogeneity and geologic age. At shorter length scales and at depths below 200 km, however, this relationship breaks down. The observed magnitude and direction of maximum anisotropy do not, in general, appear to be correlated to surface geology. The pattern of anisotropy appears to be rather complex in the upper 150 km, whereas a smoother pattern of fast axes is obtained at larger depth. If some of the deeper directions of anisotropy are aligned with the approximately N–S direction of absolute plate motion, this correspondence is not everywhere obvious, despite the fast (7 cm yr−1) northward motion of the Australian plate. More research is needed to interpret our observations in terms of continental deformation. Predictions of SKS splitting times and directions, an integrated measure of anisotropy, are poorly matched by observations of shear-wave birefringence.
1 Introduction
The Earth is heterogeneous and anisotropic on many length scales. Polarization anisotropy (modelled by transverse isotropy with a vertical symmetry axis) affects fundamental and higher-mode surface-wave propagation. Body- and surface-wave velocities are dependent on propagation azimuth. Shear-wave birefringence reflects transverse isotropy with a horizontal symmetry axis. Since anisotropy requires both the presence of anisotropic crystals and their (strain-induced) preferred orientation on seismic length scales, anisotropy measurements contain information about the mineralogy and dynamics of the mantle (see, e.g. Silver 1996; Plomerová et al. 1998; Savage 1999; Kendall 2000).
Knowledge of continental velocity heterogeneity and anisotropy may answer a number of outstanding questions. These include the very definition of the base of the continental lithosphere, the relation between the depth extent of the lithosphere and the age of the overlying crust, and the deformation of the continental lithospheric mantle and the underlying asthenosphere.
It has been argued that continental seismic anisotropy is best explained by a combination of ‘frozen’ (subcrustal) anisotropy, developed during major orogenic events, underlain by an actively deforming asthenosphere which displays anisotropy generated by the subhorizontal alignment of olivine crystals by present-day mantle flow (Silver 1996). However, considerable ambiguity exists regarding the amount and location of anisotropy inferred from body waves.
Surface-wave waveform tomography allows the delineation of lateral heterogeneities in the upper mantle with better radial resolution than body-wave traveltime tomography. A variety of methods has been developed to extract the variation of wave speed with depth from surface waves, (e.g. Nolet et al. 1986; Cara & Lévéque 1987), and a number of tomographic methods to invert for aspherical structure exist (e.g. Montagner 1986; Nolet 1990). Inversions of surface waves for heterogeneity have traditionally been restricted to fundamental modes but more recent models have incorporated higher-mode measurements, thereby enhancing the resolution down to depths of about 400 km. Inversions for both heterogeneity and azimuthal anisotropy have been made using fundamental modes (e.g. Montagner & Jobert 1988; Nishimura & Forsyth 1989) and multimode waves (Lévéque et al. 1998; Debayle & Kennett 2000a).
Local measurements of the seismic structure of the Australian upper mantle have indicated the presence of a high-velocity lithospheric lid down to ∼200 km in Northern Australia (e.g. Dey et al. 1993; Kennett et al. 1994), and to ∼300 km in Central Australia (Gaherty & Jordan 1995; Gaherty et al. 1999). Seismic anisotropy has been proposed for the upper mantle under Central Australia down to 200–300 km (Gaherty et al. 1999), and between 210 and 410 km depth under Northern Australia (Tong et al. 1994). Girardin & Farra (1998) report evidence for two-layer anisotropy in southeast Australia. The conflicting results of different shear-wave birefringence studies (Clitheroe & van der Hilst 1998; Özalaybey & Chen 1999) indicate that the anisotropy of the Australian continental mantle is complex and not yet well understood.
The skippy experiment (van der Hilst et al. 1994) has produced a unique broad-band data set of seismic waves propagating all across the Australian lithosphere, resulting in the path coverage used in the present study (Fig. 1). Waveform inversions of vertical-component fundamental and higher-mode Rayleigh-wave data have yielded detailed information on the 3-D isotropic S-wave speed variations of the continent and the adjacent oceanic areas (Zielhuis & van der Hilst 1996; Simons et al. 1999; Debayle & Kennett 2000a).
Path coverage. About 2250 vertical-component seismograms provided waveforms which were fitted to obtain the model presented in this paper. The stations belong to the skippy, kimba, agso, iris and geoscope arrays. We used earthquake locations from Engdahl et al. (1998).
Path coverage. About 2250 vertical-component seismograms provided waveforms which were fitted to obtain the model presented in this paper. The stations belong to the skippy, kimba, agso, iris and geoscope arrays. We used earthquake locations from Engdahl et al. (1998).
The azimuthal anisotropy measured by Debayle (1999) and Debayle & Kennett (2000a) changes from a regime with fast-changing directions in the upper 150 km to smoother North-South patterns below 150 km, and an analysis of Love and Rayleigh waves by Debayle & Kennett (2000b) shows polarization anisotropy extending down to about 200–250 km.
In this paper, we present a new model for the Australian upper mantle that takes into account both wave-speed heterogeneity and azimuthal anisotropy. We quantify resolution and uncertainty by means of synthetic tests, trade-off estimates, and an exact calculation of the horizontal resolution matrix. While there is a general agreement between the long-wavelength structures imaged in our present model and the one by Debayle & Kennett (2000a) (hereafter: DK2000), there are substantial differences in heterogeneity and anisotropic structure at wavelengths smaller than about 500 km even in regions for which both groups claim good resolution.
Like DK2000, in this paper we attempt to invert for 3-D structure starting from path-average models. The validity of the approximation of independent propagation of the individual modes of the seismogram has been explored theoretically (Kennett & Nolet 1990; Kennett 1995; Meier & Malischewsky 2000) and been found to be subject to certain frequency limitations, to which we adhere. As the frequency bounds for uncoupled mode propagation presume a scale length of lateral heterogeneity, DK2000 conclude one should be wary of comparing differences in our and their model on lateral length scales shorter than 300 km. DK2000 report noticeable differences from great-circle propagation through their model at frequencies >25 mHz. In the present paper, the higher-mode window of the waveforms is fit up to 50 mHz where possible. Ideally, the capability of our models to explain the observed waveforms needs to be tested with 3-D wave propagation techniques (e.g. Komatitsch & Tromp 1999).
In addition to differences in the waveform inversion technique, the tomographic inversion method used to obtain the model presented in this paper differs from the one used by DK2000. They use a form of the regionalized inversion method due to Montagner (1986), while we use the discretely parametrized, regularized inversion method introduced by Zielhuis & Nolet (1994). We compare these tomographic techniques in theory and in their application to data.
2 Inversion partitioning and data
Our method and the one of DK2000 both involve an inversion for aspherical earth structure that is partitioned into an initial waveform inversion that seeks to determine the average structure under the great-circle path between source and receiver and a subsequent tomographic inversion for 3-D Earth structure (which will be discussed in Section 4).
Our waveform inversion is due to Nolet et al. (1986). It produces a minimal set of independently constrained estimates of radial earth structure (Zielhuis 1992) that are used as data in the 3-D tomographic inversion. The independence of the linear constraints and the decorrelation of their errors is achieved by matrix diagonalization. This is the central distinguishing feature of Nolet's (1990) method. Furthermore, with this method the waveform analysis can be applied to separate analysis windows in the group velocity versus frequency space (see Simons et al. 1999, their Fig. 3), in order to separate fundamental from higher modes, where possible. The possibility of using two data windows prevents the highly energetic fundamental modes from dominating the waveform inversion. For their waveform inversion, DK2000 rely on the method by Cara & Lévéque (1987), which linearizes the waveform inversion by using the envelope of modal cross-correlation functions and fits the fundamental mode and the overtones simultaneously (Debayle 1999).
Synthetic experiment. (a) Distribution of 450 events and stations randomly distributed on a circle around Australia. (b) Input velocity heterogeneity, in per cent. (c) Input azimuthal anisotropy, direction of maximum velocity.
Synthetic experiment. (a) Distribution of 450 events and stations randomly distributed on a circle around Australia. (b) Input velocity heterogeneity, in per cent. (c) Input azimuthal anisotropy, direction of maximum velocity.
Both approaches rely on a path-average approximation and neglect mode coupling. As a result, some effects of complex source mechanisms and mantle heterogeneity are not adequately modelled. However, when applied within the frequency bounds derived for both fundamental and higher modes (Kennett & Nolet 1990; Kennett 1995) and for a large number of crossing paths, the path-average method leads to an acceptable approximation of the 3-D upper-mantle heterogeneity (Marquering & Snieder 1996).
For the present study we have expanded the data set of 1600 paths used by Simons et al. (1999) with additional data from the skippy project and from permanent observations operated by iris, geoscope and agso, and we applied the waveform inversion to a total of ∼2250 paths. In selecting the events we have avoided regions with complex source and mantle structure by limiting the data to events that originate mostly along or within the boundaries of the Indo-Australian plate (see Fig. 1).
Where data quality allows we fit the fundamental mode between 5–25 mHz (40–200 s period) or a narrower frequency band within this range, while the higher modes are fitted between 8–50 mHz (20–125 s) or less. With these frequency bounds the fundamental modes can provide sensitivity to structure down to at least 400 km, and the higher modes yield additional detail. For more details about the application to skippy data, see Zielhuis & van der Hilst (1996).
At high frequencies, Rayleigh modes are sensitive to variations in crustal thickness and velocity. We use continental crustal thickness determinations from broad-band receiver functions at each of the skippy and permanent network stations (Shibutani et al. 1996; Clitheroe et al. 2000a,b). This is arguably the most accurate crustal model to date of the Australian continent.
3 Azimuthal anisotropy
In the tomographic inversion we will take as a starting point the particular radial velocity averages that result from the waveform inversion (see Section 2), and solve for heterogeneity (βV) and the anisotropic parameters (A1=Gc/(2ρβV) and A2=Gs/(2ρβV)). We calculate the anisotropic magnitude as (G2c+G2s)1/2 and the direction of maximum velocity as tan−1(A2/A1)/2.
The ubiquitous use of the parametrization represented by eqs (1) and (6) is based on the dominance of the ∂cR/∂L partial for fundamental-mode surface-wave propagation (Montagner & Nataf 1986). Since ∂cR/∂L is also dominant for higher modes, both for isotropic and realistically anisotropic media (to verify this, we computed sensitivity kernels for the model by Gaherty & Jordan 1995), a parametrization in terms of βV, Gc and Gs is adequate when higher modes are included.
To fit seismograms, we used the ∂cR/∂L partial derivative with an isotropic reference model. The resulting 1-D path model of βV is an adequate representation of the Earth when no further azimuthally anisotropic information is available, because the sensitivity kernels of Rayleigh wave phase velocities to perturbations in S-velocity in isotropic and transversely isotropic earth models are very similar (Anderson & Dziewonski 1982), although a small error is probably made (Lévéque & Cara 1985).
4 The tomographic problem
4.1 Construction of a linear matrix problem
4.2 Model parametrization
The model presented in this paper is based on a 2°× 2° equal-area surface parametrization of the kind shown in Figs 2(a) and (b). An inversion on a 3°× 3° grid showed essentially the same features, but an inversion for heterogeneity and anisotropy on a 1°× 1° grid proved less stable. To further test the robustness of the structures retrieved, we performed inversions on a resolution-dependent grid (Fig. 2c). In order to determine the extent to which the Australian continental structure at varying depths bears a relationship to the geology observed at the surface, we also performed inversions on a pure-path (tectonic) grid (Fig. 2d).
(a) Path density and (b) azimuthal coverage for the paths shown in Fig. 1, calculated on a 2°× 2° equal-area grid. As a measure of azimuthal coverage, we took the number of pairs of rays in the cell coming into two different quadrants. (c) Design of a grid adapted to expected model resolution. Path density and azimuthal coverage calculated on a 1°× 1° square grid were normalized and added to provide a measure of data quality. The cell size was adapted to allow a finer parametrization in areas where the resolution is expected to be better. (d) Tectonic subdivision. Based on an underlying equal-area grid, similarly shaded regions indicate approximate tectonic domains which are used as a single parameter in the pure-path inversion. Lines indicate the location of major sedimentary basins. For the definition of tectonic units, see Simons et al. (1999).
(a) Path density and (b) azimuthal coverage for the paths shown in Fig. 1, calculated on a 2°× 2° equal-area grid. As a measure of azimuthal coverage, we took the number of pairs of rays in the cell coming into two different quadrants. (c) Design of a grid adapted to expected model resolution. Path density and azimuthal coverage calculated on a 1°× 1° square grid were normalized and added to provide a measure of data quality. The cell size was adapted to allow a finer parametrization in areas where the resolution is expected to be better. (d) Tectonic subdivision. Based on an underlying equal-area grid, similarly shaded regions indicate approximate tectonic domains which are used as a single parameter in the pure-path inversion. Lines indicate the location of major sedimentary basins. For the definition of tectonic units, see Simons et al. (1999).
Due to the uneven sampling density the resolution is not uniform across the model domain, but the grid spacing can be adapted to the locally expected resolution in order to reduce the number of free parameters. Our approach to adaptive grids is akin to that of Abers & Roecker (1991). We used a fine (1°× 1°) Cartesian grid upon which we define a coordinate transform that attempts to equalize the sampling density by clustering cells together. Starting from a maximum cell size of 8°× 8° we allow a further subdivision, into cells sized 4°× 4°, 2°× 2° and 1°× 1°, as long as this brings the resulting mean of the sampling density over these new cells closer to the global mean. As a measure of resolution we use the sum of path density (see Fig. 2a) and azimuthal variability (Fig. 2b) in the surface cells. Although this is not necessarily an optimal a priori estimate of the spatial resolution of the model the similarity of the features of this measure with an exactly computed resolution matrix (see Section 5.2) is sufficient to warrant their use in the construction of an irregular grid. An example of an irregular grid that results from the clustering of the sampling density is plotted in Fig. 2(c).
Pure-path methods are used to invert for the average velocity and anisotropy structure of predefined tectonic regions (Fig. 2d) (Nishimura & Forsyth 1985). Similarly to the construction of a resolution-dependent irregular grid we achieve this tectonic parametrization by a coordinate transform, such that all surface cells belonging to a tectonic domain are treated as a single parameter in the inversion.
4.3 Regularization or a priori information?
The discrete surface parametrization of our method leads to classic regularization approaches of systems of linear equations (Section 4.3.1). On the other hand, the continuous regionalization approach adopted by DK2000 draws solutions from a class of functions defined by their a priori covariance structure (Section 4.3.2).
4.3.1 Discrete regularization
The optimal weighting matrix A−1 is given by the inverse of the data covariance matrix Cd−1, so that the solution has minimum variance. A covariance matrix, Cd is positive definite (hence invertible) and symmetric; in our implementation it is also diagonal because the individual parameters resulting from the waveform fitting procedure are independent by construction (Nolet 1990; Zielhuis & Nolet 1994).
If we choose B−1=I, the identity matrix, we are imposing a minimum norm constraint on the solution by biasing the solution to zero in a least-squares sense. This is appropriate if m represents perturbations to a reference model which is sufficiently close to the real solution so that linearity is preserved. Most inversion algorithms for the generalized inverse minimize the model norm implicitly (van der Sluis & van der Vorst 1987). In particular, the LSQR algorithm provides an intrinsic damping of the norm m·m when the iterations are truncated at some small finite value (Nolet 1985; van der Sluis & van der Vorst 1987).
The effect of this “roughening” operator is to constrain the solution to be locally smooth (Nolet 1987). Eq. (11) is the Toeplitz matrix of a high-pass filter in model space (Strang & Nguyen 1997). Hence, we are imposing that the high-passed model vector be minimal, which, after inversion, results in a smooth model.
4.3.2 Continuous regionalization
4.3.3 Equivalence of regularization and Bayesian approaches
Regularization of a discretely parametrized model and stochastic extensions using an a priori model covariance structure have been presented as distinct alternatives (see, e.g. Scales & Snieder 1997) but the approaches are equivalent in many ways (e.g. Ho-Liu et al. 1989).
With this choice of covariance function one obtains smooth, infinitely differentiable functions. Eq. (20) corresponds to a model norm weighted by an infinite number of finite difference matrices, starting with the zeroth difference, i.e. the traditional model norm. Neglecting the higher-order smoothing terms makes this approach equivalent to eq. (19).
To illustrate this equivalence we conducted a synthetic experiment to compare the discrete regularization approach with the continuous regionalization for one particular depth. For a fabricated path coverage (Fig. 3a), (noiseless) data were synthesized from a model of shear-wave speed heterogeneity (Fig. 3b), azimuthal anisotropy (Fig. 3c), or both. Fig. 4 shows the inversion results obtained with the discrete method (left column) and the continuous method (right column). Figs 4(a) and (b) show the result for an input model with only heterogeneity (i.e. Fig. 3b). Figs 4(c) and (d) show the solution for a model with only anisotropy (i.e. Fig. 3c), and Figs 4(e) and (f) give the results for a model which has both heterogeneity and anisotropy. This test demonstrates that a set of inversion parameters (correlation lengths, a priori standard deviations and regularization parameters) can be found with which both methods yield very similar results. The exact values are unimportant given the synthetic nature of this example. Differences in amplitude recovery are due to the different implicit weighting of the data constraints versus the a priori information in both methods. The differences in smoothness are due to the chosen horizontal correlation length (chosen to maximize variance reduction; in this example we chose a correlation length of 175 km for both the velocity heterogeneity and the anisotropy) and the implicit correlation length provided by the damping.
Comparison of the discrete iterative approach (left column) with the continuous regionalization inversion (right column). (a, b) Isotropic inversion result with the path distribution of Fig. 3(a) and input heterogeneity of Fig. 3(b). (c, d) Anisotropy-only inversion with input anisotropy of Fig. 3(c). (d, e) Joint inversion for velocity heterogeneity and anisotropy with input anomalies as in Figs 3(b) and (c).
Comparison of the discrete iterative approach (left column) with the continuous regionalization inversion (right column). (a, b) Isotropic inversion result with the path distribution of Fig. 3(a) and input heterogeneity of Fig. 3(b). (c, d) Anisotropy-only inversion with input anisotropy of Fig. 3(c). (d, e) Joint inversion for velocity heterogeneity and anisotropy with input anomalies as in Figs 3(b) and (c).
4.4 Differences between Method I and Method II
While conceptually similar there are, however, theoretical and practical differences that set our method (Method I) apart from the method used by DK2000 (Method II).
In the previous section we have shown the theoretical equivalence of the Bayesian inversion approach with the regularization method, and that inversion parameters can be found so that both methods produce similar results (Figs 3 and 4). One should realize, however, that we could find such parameters only because the input model was known. In inversions for unknown earth structure this tuning is not possible. The challenge is to find a priori horizontal correlation lengths for both the velocity heterogeneity and the anisotropic parameters, and a priori standard deviations of the anomalies (σβ for δβV and σA for A1 and A2). The correlation length affects the smoothness of the resulting model and the a priori uncertainties regulate the amplitude. Moreover, the ratio of σβ and σA controls the trade-off between these quantities. The choice of σ does not affect the location but only the amplitude of the anomalies. Various authors using this technique have made different choices. Lévéque et al. (1998), for example, choose σβ= 100 m s−1 and σA= 5 m s−1. Griot et al. (1998) have σβ= 200 m s−1 and σA= 200 m s−1, and Silveira et al. (1998) take σβ= 200 m s−1 and σA= 130 m s−1. DK2000 take σβ= 50 m s−1 and σA= 3 m s−1.
In contrast, in Method I we incorporate the sensitivity to heterogeneity and anisotropy with an equal weighting. The relative amplitudes of anisotropy and heterogeneity are controlled by the weights of the damping terms of their norm and not well constrained by the data. In our model, the norm damping for the anisotropic parameters was twice as large as for the velocity heterogeneity. We evaluate the amount of trade-off a posteriori in various experiments (see Section 5.5).
Another difference is the treatment of the third dimension in the 3-D inversion. In Method I the third dimension is an integral part of the minimization problem, entering the equations through the depth-dependent linear constraints on the velocity and anisotropic structure (see Zielhuis & Nolet 1994). In this method, the path-averaged constraints entering the tomographic inversion are depth-dependent functions representing the particular depth-average of structure that can be resolved independently from all other constraints (Nolet 1990; Zielhuis 1992). With the continuous regularization technique of Montagner (1986) such a treatment of the depth-dependent resolvability of Earth structure can be achieved by retaining the full 3-D covariance matrices C(r, r′) in the development. However, DK2000 invert the shear-wave speed path-averages layer by layer. The data used in their tomographic inversion are given by the shear-wave speed profiles obtained from the waveform inversion, sampled with an interval of 25 km, and individual depths are treated separately. The results are then stacked into a 3-D model. The variation of the resolution matrix with depth is given by the errors on the shear-wave speed at the particular depths. However, judging from the kernel functions that represent how surface waves sample the Earth, it is not possible to determine shear-wave speeds at such finely spaced depth independently from one another. As a consequence, with this method, it is impossible to obtain direct estimates of the variation of lateral resolution with depth and the extent of radial interdependence in the 3-D model.
Finally, in Method I, we use the inversion algorithm LSQR (Paige & Saunders 1982), which speeds up the computations by several orders of magnitude, but no full resolution matrix can be computed directly. In contrast, Method II performs the inversion directly, which is computationally much more expensive but enables the explicit computation of the resolution matrix, a capability we will use in the next section.
5 Resolution and robustness analysis
5.1 Data coverage
The best sampled regions are the easternmost two-thirds of the Australian continent (Fig. 1), with a maximum path density on the eastern rim of the continent (Fig. 2). Towards the west, the path density degrades. Fig. 2(b) depicts the number of pairs of rays that come into the cell in a different azimuthal quadrant. This score is more uniform across the continent, but it also diminishes in the western half. The northeastern portion of the continent has a maximum of about 200 such pairs of crossing rays in one cell. In Fig. 2(c) we combined both resolution criteria into a single value ranging from 0 to 1.
5.2 Exact calculation of the resolution matrix
Using the direct Bayesian inversion we calculate the exact resolution matrix both for the actual (Fig. 1) and the synthesized uniform (Fig. 3a) data coverage and compare this to the qualitative resolution estimates from checkerboard tests. We input the path-average velocity anomaly at a single depth into the Bayesian inversion code, with a correlation length of 175 km for both velocity and anisotropy, and a priori model uncertainties of 120 m s−1 for δβV, A1 and A2 (see eq. 7). The value of these choices only affect the amplitude of the result; equating 
and σA corresponds to an equal weighting of velocity heterogeneity and anisotropy.
The results for the synthetic uniform path coverage of Fig. 3(a) are plotted in Fig. 5(a) (resolution of δβV), Fig. 5(b) (resolution of A1), and for the actual coverage of Fig. 1 in Figs 5(c) and (d). The resolution for A2 is virtually identical to that of A1 and is therefore omitted. The comparison of the synthetic resolution with the resolution from the actual data gives us a feel for the values we can expect where the model is well resolved. As anticipated, the solution in the eastern two-thirds of the continent is markedly better than in the west.
Resolution tests. Exact horizontal resolution computation for the synthetic path density experiment (a–b) and the real data coverage (c–d). Plotted is the diagonal value of R(r, r′) for the isotropic velocity solution δβV (a, c), and the cosine term A1 (b, d). The resolution for the sine term A2 is virtually indistinguishable from that of the cosine term.
Resolution tests. Exact horizontal resolution computation for the synthetic path density experiment (a–b) and the real data coverage (c–d). Plotted is the diagonal value of R(r, r′) for the isotropic velocity solution δβV (a, c), and the cosine term A1 (b, d). The resolution for the sine term A2 is virtually indistinguishable from that of the cosine term.
5.3 Checkerboard tests
A popular way of assessing the resolution in tomographic models is to calculate the recovery of a synthetic (e.g. checkerboard) input pattern. By themselves checkerboard tests are difficult to interpret. It is often hard to guess the exact pattern of the anomalies without knowledge of the input function, which strongly degrades their diagnostic value. Also, the damping parameters can be chosen to reproduce the input pattern optimally, a luxury which is not available when choosing the damping needed to model the actual data for unknown Earth structure. These and other, intrinsic, limitations of checkerboard tests have been pointed out by several authors (e.g. Lévéque et al. 1993). Hence, we will use the results of checkerboard tests in a qualitative way only, and compare them to the exact resolution computed in Section 5.2.
Since the sensitivity matrix maps the input anomalies onto the synthetic data according to the independently constrained depth-dependent kernel functions, the solution at one particular depth is related to the solution at adjacent depths in Method I. Therefore, we can use checkerboard tests to get a sense of the radial resolution and vertical smearing of the tomographic model (Fig. 6) as well as the horizontal resolution and its depth variation (Fig. 7).
Checkerboard test for vertical smearing. An input pattern of velocity heterogeneity at 80 km (a) superimposed on a pattern of azimuthal anisotropy (b; shaded for clarity) offset by half a cell was used to generate synthetic data. The inversion for both wave speed anomalies and anisotropy yields a recovery at 80 km given in (c) and (d). Due to the radial averaging properties of the surface-wave kernels and the smearing in the inversion, some of the input structure at 80 km is also mapped into the next layer at 140 km (e and f). Panels a, c, and e are plotted on the same color scale.
Checkerboard test for vertical smearing. An input pattern of velocity heterogeneity at 80 km (a) superimposed on a pattern of azimuthal anisotropy (b; shaded for clarity) offset by half a cell was used to generate synthetic data. The inversion for both wave speed anomalies and anisotropy yields a recovery at 80 km given in (c) and (d). Due to the radial averaging properties of the surface-wave kernels and the smearing in the inversion, some of the input structure at 80 km is also mapped into the next layer at 140 km (e and f). Panels a, c, and e are plotted on the same color scale.
Checkerboard test for depth-varying horizontal resolution. (a, b) Percentage of input anomaly recovered of an isotropic checkerboard pattern at 80 and 200 km depth. (c, d) Difference in angle between input and output anisotropy in a joint inversion for a pattern of heterogeneity and anisotropy as in Figs 6(a)–(b) at 80 and 200 km depth. Bars scale with the magnitude of the difference such that excellent recovery is depicted by a dot.
Checkerboard test for depth-varying horizontal resolution. (a, b) Percentage of input anomaly recovered of an isotropic checkerboard pattern at 80 and 200 km depth. (c, d) Difference in angle between input and output anisotropy in a joint inversion for a pattern of heterogeneity and anisotropy as in Figs 6(a)–(b) at 80 and 200 km depth. Bars scale with the magnitude of the difference such that excellent recovery is depicted by a dot.
We have created two staggered checkerboard patterns, one for heterogeneity and one for azimuthal anisotropy (Figs 6a and b). Synthetic data were generated for a model with this input pattern introduced at 80 km depth. The inversion result is shown at this same depth (80 km: Figs 6c and d) as well as at the next level of the solution (140 km: Figs 6e and f). Except for the diminished amplitudes the velocity input pattern is reasonably well recovered, including the sharp gradients between the cells. This and other tests suggest that large wave speed contrasts occurring over as little as 250 km can be resolved by our data. On the other hand this test does not tell us the amplitude of the minimum wave speed variation that can be detected at that same wavelength.
The N–S and E–W pattern of azimuthal anisotropy is fairly well recovered in the central and eastern part of the continent but the angular variability is at least 20° (Fig. 6d). Indeed, in large parts of our study region it is not obvious from Fig. 6(d) alone that the input pattern was a regular checkerboard pattern, which hampers the interpretation of the results of such tests significantly.
Radial smearing affects our results, but a comparison of the amplitudes of the solutions at both depths reveals that the level of contamination is small.
In Fig. 7 we illustrate the variation of horizontal resolution with depth. We have plotted the output-to-input ratio (in per cent) of an isotropic checkerboard test with patches of about 550 km each (5°× 5°, performed at 80 (Fig. 7a) and 200 km (Fig. 7b)). At 80 km, we observe that up to 80 per cent of the input anomaly is recovered in the best resolved regions; at 200 km this number has dropped to about 60 per cent, but in most of the study region the amplitude recovery is less than 50 per cent. Figs 7(c) and 7(d) plot the difference in angle between the synthetic input pattern of Fig. 6(b) in the presence of the velocity heterogeneity of Fig. 6(a) and the recovery of the anisotropy from a joint inversion for heterogeneity and anisotropy, for an input pattern and model retrieval at 80 (also shown in Fig. 6d) and 200 km depth (also shown in Fig. 6e).
We conclude that, except perhaps in western Australia, our final model will give a fair representation of the velocity and anisotropy structure of the upper mantle under Australia.
5.4 Model robustness
In addition to the parametrization strategies described in Section 4.2 on the grids shown in Fig. 2, we evaluate the robustness of the models by an extensive series of tests. We performed inversions with up to 40 per cent of the data randomly removed. We also calculated the mean model obtained from 500 inversions with 5 per cent of the data randomly removed, as well as the mean of 500 inversions performed on a data set to which noise was added (with signal-to-noise ratios up to 15 dB). The average models obtained in this way were visually indistinguishable from our final model. For this class of mixed-determined inversion problems it is not possible to use inversion experiments like these to interpret the variance of the model parameters quantitatively everywhere. Areas with low resolution are always damped towards the reference model, and this perfectly known state of total uncertainty has therefore almost zero variance in such inversions. However, for the well resolved areas in the east, the standard deviation obtained is everywhere less than 0.8 per cent, which is much smaller than the magnitude of the anomalies we obtain in the final model. This bound holds for all depths of our final model.
5.5 Trade-off of heterogeneity and anisotropy
Several previous studies (e.g. Vinnik et al. 1992; Gaherty & Jordan 1995; Clitheroe & van der Hilst 1998; Girardin & Farra 1998; Debayle & Kennett 2000a,b) have suggested that the upper mantle beneath Australia is anisotropic. However, we need to show that the inclusion of the azimuthal anisotropy into the inversion problem is required by our data.
The increase in variance reduction of isotropic-only (82.6 per cent in a separate inversion with all parameters unchanged) versus azimuthally anisotropic models (86.1 per cent for our model) is an indication of the improvement of data fit when anisotropy is included, but its statistical significance is hard to assess. Methods to evaluate the improvement of data fit with an increasing number of parameters, such as Fischer tests, are only applicable to overdetermined problems (Menke 1989) and can thus not be used here. One indication of improvement of data fit is by comparing the increase in variance reduction from an isotropic model (δβV) parametrized in 7059 surface cells of a certain area with an anisotropic inversion (δβV, A1, A2) of 2353 cells of an area three times as large. The number of inversion parameters is the same in either case, but the variance reduction increases from 83.5 to 86.1 per cent. Thus, the inclusion of anisotropy results in an increase of variance reduction even when the number of free parameters is left unchanged.
We have studied the incremental improvement in variance reduction in our final model due to the inclusion of progressively deeper layers of heterogeneity, anisotropy, or both. For the velocity heterogeneity alone, the depth intervals between 30 and 400 km depth contribute most to explain the variance in the data, reaching a plateau of 76.9 per cent for the model between 0–400 km and a maximum of 77.5 per cent for the entire depth range (0–670 km). Anisotropy alone between 0–400 km only explains 20.1 per cent of the variance in the data, terminating at 20.5 per cent for the entire model between 0–670 km. The improvement in variance reduction by the addition of successive layers of anisotropy to a model with heterogeneity between 0–670 km is from 77.5 per cent to 85.9 per cent at 400 km, terminating at 86.1 per cent for the entire model.
We explore the significance of the azimuthal anisotropy with another synthetic experiment designed to test to which extent uneven data coverage and heterogeneity can combine to produce an anisotropic bias. We represent our full solution by (δβV, A1, A2). In one test we inverted data synthesized from the heterogeneity part only, (i.e. (δβV, 0, 0)) for both heterogeneity and anisotropy. From the result we can then infer how much anisotropy is produced by uneven data coverage and heterogeneity. In a second test we invert synthetic data from (0, A1, A2) for heterogeneity and anisotropy to assess the projection of artificial anisotropy into heterogeneity.
Fig. 8 summarizes the test results for a depth of 140 km. Figs 8(a) and 8(b) show the inversion of models which should have zero residual anisotropy and zero residual velocity heterogeneity, respectively. The magnitude of the anisotropy in the actulal model plotted in Fig. 8(c), can be compared to the magnitude of the “erroneous” anisotropy (Fig. 8d). For all depths considered, the average anisotropic bias is less than 20 per cent of the magnitude required to satisfy the data in the actual model.
Trade-off of heterogeneity with anisotropy at 140 km. (a) Inversion for heterogeneity and anisotropy of synthetic data calculated from the isotropic part of our model. With perfect resolution and control of the trade-off, the anisotropy retrieved should be zero everywhere. (b) Inversion of synthetic data calculated from the anisotropic part of the model for heterogeneity and anisotropy. In this case perfect resolution should have resulted in zero heterogeneity. (c) Magnitude of anisotropy (in GPa) of the actual model (see Fig. 9). (d) Magnitude of the erroneous anisotropy of Fig. 8(a). Contour lines plotted at 1 GPa interval for both panels.
Trade-off of heterogeneity with anisotropy at 140 km. (a) Inversion for heterogeneity and anisotropy of synthetic data calculated from the isotropic part of our model. With perfect resolution and control of the trade-off, the anisotropy retrieved should be zero everywhere. (b) Inversion of synthetic data calculated from the anisotropic part of the model for heterogeneity and anisotropy. In this case perfect resolution should have resulted in zero heterogeneity. (c) Magnitude of anisotropy (in GPa) of the actual model (see Fig. 9). (d) Magnitude of the erroneous anisotropy of Fig. 8(a). Contour lines plotted at 1 GPa interval for both panels.
In regions of uneven data coverage some artificial anisotropy is clearly produced. The direction of the artificially induced anisotropy is especially sensitive to uneven sampling, but, in general, the magnitude of the bias is small. In Western Australia, however, much of the anisotropy present in the actual model may well consist of bias due to the scarcity of crossing paths in that region. This makes it impossible to distinguish the fast velocity anomaly, which is clearly detectable in the waveforms, from a possible alignment of the fast axis of anisotropy.
6 Results
In Figs 9 and 10 we present models for heterogeneity and azimuthal anisotropy of Australia, using Method I. The inversions were performed on the regular (equal-area) grid shown in Figs 2(a) and (b).
Inversion results at different depths. Percentages from reference model with background shear-wave velocities of 4500 m s−1 at 80 and 140 km depth; 4513.5 m s−1 at 200 km; 4635 m s−1 at 300 km.
Inversion results at different depths. Percentages from reference model with background shear-wave velocities of 4500 m s−1 at 80 and 140 km depth; 4513.5 m s−1 at 200 km; 4635 m s−1 at 300 km.
Cross-sections at various latitudes. Fifteen contour lines are plotted between −6 and 6 per cent.
Cross-sections at various latitudes. Fifteen contour lines are plotted between −6 and 6 per cent.
The azimuthally isotropic component δβV shows many of the characteristics that were present in our earlier models (Zielhuis & van der Hilst 1996; Simons et al. 1999). Low wave speed anomalies mark the eastern Phanerozoic rim of the continent. They are associated with thermal anomalies, recent volcanism and high attenuation (Mitchell et al. 1998). Wave speeds in the central and western region are high, but the correlation between the high velocities and the age of the overlying crust is scale-dependent (Simons et al. 1999).
In most regions the anisotropy suggested by the data exceeds the anisotropic bias (see also Fig. 8), and the lateral variability in directions is much larger than the nominal uncertainty of 20°. There is no obvious relationship between wave speeds and anisotropy, neither in magnitude nor in direction, nor is there a clear correlation between the anisotropy and the large-scale geological structure (Fig. 11).
Average heterogeneity and anisotropy in function of crustal age. A coarse regionalization of the Australian continent consists of a predominantly Archaean domain (A), the central Proterozoic region (Pr) and domains of Phanerozoic age (Ph1 and Ph2). While wave speed anomalies become progressively more positive with increasing crustal age (see also Simons et al. 1999), no defining relationship is found in the magnitude of the anisotropy.
Average heterogeneity and anisotropy in function of crustal age. A coarse regionalization of the Australian continent consists of a predominantly Archaean domain (A), the central Proterozoic region (Pr) and domains of Phanerozoic age (Ph1 and Ph2). While wave speed anomalies become progressively more positive with increasing crustal age (see also Simons et al. 1999), no defining relationship is found in the magnitude of the anisotropy.
The cross-sections shown in Fig. 10 suggest a depth to the base of the continental lithosphere of up to 250–300 km (see also Simons et al. 1999; Fischer & van der Hilst 1999), and Fig. 11 shows that the wave speed difference between the large-scale tectonic provinces is confined to the shallowest 300 km. The central-Australian region is underlain by the thickest high-velocity lid, which may be an effect of preferential preservation (Simons et al. 1999). Even with the uncertainty estimates given in Section 5.5, the magnitude of the anisotropy does not fall off with depth as rapidly as the magnitude of the velocity heterogeneities. The average anisotropic magnitudes reach a maximum around 100 km (Fig. 11) and diminish at larger depth, in accordance with observations elsewhere (Babuška et al. 1998). For comparison, the magnitudes of the anisotropy from laboratory estimates of upper-mantle materials range from 2 to 7 GPa (Peselnick & Nicolas 1978; Estey & Douglas 1986).
Although not all individual directions of anisotropy agree with the model by DK2000, we concur with these authors that the pattern of anisotropy changes with increasing depth. The solutions at 80 and 140 km depth show a complex pattern of azimuthal anisotropy (short correlation wavelengths), whereas in the layers at 200 and 300 km a simpler, smoother pattern (long correlation wavelengths) is observed (Fig. 9). The character of the anisotropy and the length scale over which it varies seem to change rapidly across the 150 and 200 km depth interval.
The velocity heterogeneity revealed in Fig. 9 was reproduced satisfactorily with the irregular parametrization (Fig. 2c), albeit with somewhat reduced amplitudes. In the western regions, where the path coverage is demonstrably worse, our fine equal-area grid over-parametrizes the structure but this does not introduce noticeable artifacts. Within the uncertainty of 20°, the different parametrizations yield similar directions of fast axes of anisotropy. The pure-path average solutions (on the grid shown in Fig. 2d) at 80, 140 and 200 km were very similar to the averages of the wave speeds imaged on the equal-area grid. Differences appear at larger depth, suggesting again hat below 200 km most correspondence between wave speed structure and surface geology is lost. The central and northern region of Australia have stable directions of maximum anisotropy, both under the equal-area and the tectonic parametrization, but the magnitudes are much reduced in the latter. In the eastern portion of the continent a tectonic parametrization reveals little anisotropy, at any depth. This could indicate that the anisotropy is small, i.e. at or below the level of the anisotropic bias, or that through much of the depth range the anisotropy varies on length scales smaller than the main tectonic units of the Australian continent.
7 Discussion
7.1 Thickness, age-relations and deformation of the Australian lithosphere
Many global models of Rayleigh wave phase velocities (e.g., Trampert & Woodhouse 1996; Ekström et al. 1997) have measured faster than average wave speeds for the Archaean and Proterozoic portions of the Australian continent (which correspond roughly to the western two-thirds of the continent). This pattern is corroborated by our regional model, but significant variations in heterogeneity and anisotropy exist at smaller scales than is apparent from global models. In particular, the depth extent of the lithosphere is determined in much more detail than possible with global data sets. Some global models suggest that the high wave speeds associated with ancient continents extend to depths exceeding 400 km (e.g. Su et al. 1994; Masters et al. 1996). Our detailed 3-D model verities instead the local measurements of lithospheric thickness of ∼200 km quoted for Northern Australia (Dey et al. 1993; Kennett et al. 1994; Tong et al. 1994), and ∼300 km under Central Australia (Gaherty & Jordan 1995; Gaherty et al. 1999). A ∼200 km thick lithosphere in North Australia is compatible with our observations (see, e.g. Fig. 10a), but clearly the 3-D structure of the lithosphere–asthenosphere boundary is complex, and the thickness of the lithosphere probably reaches at least 300 km in certain parts of Central Australia (Fig. 10b). Eastern Australia, which is the most recently tectonized region of the continent, is characterized by a pronounced low velocity zone between 100–200 km depth underlain by fast anomalies at a depth of about 300 km (see Fig. 10b). The origin of the latter is uncertain and it may not be related to the overlying lithosphere. We conclude that the thickness of the Australian lithosphere is about 225 ± 50 km. This is in agreement with recent observations for the Canadian Shield (Frederiksen et al. 2001) and South-Africa (James et al. 2001). Similar lithospheric thicknesses are obtained for other cratonic regions (for a review, see Artemieva & Mooney 2001). Elsewhere, Simons & van der Hilst (2002) report the full variation of lithospheric thickness based on contours of the velocity anomalies in our present model in more detail.
The wave speed heterogeneity of the Australian continent seems to correlate with crustal evolution and lithospheric age, at least to depths of about 200 km (Fig. 11) (see also Simons et al. 1999). At depths greater than 250 km there is no unambiguous relationship between crustal age and wave speed. The anisotropy does not seem to change on the length scale of the major geological domains at any depth, although its magnitude remains higher at larger depth beneath the regions affected by the most recent tectonic activity (the easternmost parts of the continent). In contrast, Babuška et al. (1998) argued that there are systematic variations in radial anisotropy which are related to the age of continental provinces in North America and Eurasia.
Global and regional studies (e.g. Montagner & Tanimoto 1991; Lévéque et al. 1998) have been interpreted to show evidence of a shallow regime of rapidly varying anisotropy underlain by a region with smoother patterns of anisotropy. DK2000 observe this change at 150 km depth. There is uncertainty on the magnitude and directions of the azimuthal anisotropy, but qualitatively, the wave speed and anisotropy maps presented in Fig. 9 argue for a change in regime between the first two depths (80 and 140) and the deeper parts of the solution (200 and 300 km). The first regime is characterized by a complex pattern of strong anisotropy, whereas the second presents smoother variations of weaker anisotropy and a tendency towards N–S fast axes increasing with depth.
Gaherty et al. (1999) and Debayle & Kennett (2000b) found that polarization anisotropy in the Australian continent is concentrated above 250 km. Our Rayleigh wave model does not contain information about polarization anisotropy. However, most of the anisotropic signal is present in the depth range above 200 km depth (Fig. 11), and the improvement in variance reduction due to the inclusion of azimuthal anisotropy is largely explained by the layers above 300 km depth. While we believe we see evidence of azimuthal anisotropy until a depth of about 300 km, its magnitude is not incompatible with the statement by Gaherty et al. (1999) that polarization anisotropy below 250 km depth must be small (<1 per cent). We interpret our model to show evidence for fossil anisotropy above 200 km and deeper anisotropy which can be reconciled with the overall direction of absolute plate motion. The precise interpretation of this model in terms of continental deformation will be the subject of forthcoming research.
7.2 Comparison with results from shear-wave splitting
From the analysis of refracted S-waves bottoming under North Australia, Tong et al. (1994) obtain shear-wave splitting times exceeding 4 s in certain places. Vinnik et al. (1989) and Barruol & Hoffman (1999) report no perceptible splitting for SKS arrivals at station CAN in southeastern Australia, but Girardin & Farra (1998) interpret P-to-S converted phases to be indicative of a two-layer anisotropic system under that station with nearly perpendicular fast axes distributed over comparable depth intervals so the contributions to splitting cancel out. Clitheroe & van der Hilst (1998) present evidence for complex, frequency-dependent splitting from SKS and SKKS phases, with split times generally less than 1 s in the Precambrian central part of Australia. Their fast axis directions do not align with the direction of absolute plate motion, which is roughly SSW-NNE, but rather suggest alignment with crustal fabric or a complex lattice preferred orientation of olivine in the subcrustal lithosphere. However, Özalaybey & Chen (1999) find no convincing evidence for splitting at all. They argue that the analysis of a small amount of data over a limited azimuthal range and complex multipathing effects due to heterogeneities at depth may contribute to erroneous splitting measurements, especially at the frequencies higher than 0.3 Hz considered by Clitheroe & van der Hilst (1998).
Applying eqs (22) and (23) to our model (depicted in Fig. 9) yields considerable splitting times throughout the continent (Fig. 12a). The results of the surface-wave inversion presented here suggest higher levels of azimuthal anisotropy than inferred from SKS by either Clitheroe & van der Hilst (1998) or Özalaybey & Chen (1999). For comparison, we have plotted a recent data set of splitting times and directions in Fig. 12(b) (from Iidaka, in prep.).
(a) SKS splitting directions predicted from the solution given in Fig. 9. Station CAN is represented by the white triangle. (b) Shear-wave splitting observations (from Iidaka, in prep.).
(a) SKS splitting directions predicted from the solution given in Fig. 9. Station CAN is represented by the white triangle. (b) Shear-wave splitting observations (from Iidaka, in prep.).
In addition to calculating SKS times due the entire depth range of anisotropy, we have investigated the incremental effect on the splitting times from successively deeper levels in our model. The transverse isotropy (with horizontal symmetry axis) between 140 and 400 km depth contributes most to the calculated SKS split times. This is in good agreement with Tong et al. (1994), who argue that, in Northern Australia, the most likely source region of anisotropy is the asthenospheric zone between 210–410 km.
The fast axis directions in the eastern portions of the continent follow the rough outline of the high wave speeds associated with the Australian Precambrian (Fig. 9a). Clitheroe & van der Hilst (1998) have suggested this might be related to the crystal alignment by asthenospheric flow around the irregularly shaped eastern edge of the deep continent (see also Fouch et al. 2000). This is supported by the fact that the split times originate mainly at a depth near 140 km, in a zone of anomalously low velocities and, presumably, low viscosities (Figs 10 and 11). In the center of the continent, the predicted splitting times are somewhat smaller, which may in part be due to a change from E–W in the shallower part to N–S in the deeper continental lithosphere (see, e.g. Silver & Savage 1994; Saltzer et al. 2000).
A mismatch between surface-wave and body-wave anisotropy has been observed on other continents as well (Saltzer et al. 2000; Freybourger et al. 2001). We should point out that the assumption of a horizontal orientation of the fast axis is not unlikely to be violated, especially in the shallow parts of the continental lithosphere. Further analysis will require taking into account the depth variation of fast axis directions, the possibility of dipping axes and frequency-dependent effects.
Our observations point to a complex geometry of the continental anisotropy at various length scales. We cannot rule out that strong anisotropy in the top 80–140 km exists if the direction of fast axes varies rapidly on a scale that is smaller than can be resolved by our data. With the magnitude of anisotropy increasing and the lateral correlation lengths decreasing towards the surface, it is conceivable that strong, rapidly changing azimuthal anisotropy in the top 140 km diminishes the effective birefringence of vertically propagating body waves without similarly influencing the long-period surface-wave data used in our study (Saltzer et al. 2000).
8 Conclusions
Manual data selection and processing of stations operated by agso, iris and geoscope and data from the portable skippy experiments have produced a data coverage that is superior to surface-wave studies on any other continent in terms of density and azimuthal distribution.
We have compared two tomographic approaches for the inversion of path-average earth models for 3-D models of S-heterogeneity and azimuthal anisotropy. Method I (Zielhuis & van der Hilst 1996; Simons et al. 1999) performs a 3-D inversion for heterogeneity and anisotropy based on a discrete parametrization at the surface with depth basis functions representing independently constrained estimates of radial earth structure. Method II, used by Debayle & Kennett (2000a), uses a continuous regionalization algorithm which performs a separate inversion for each depth layer. Discretely parametrized regularization methods and Bayesian continuous regionalization methods can lead to similar results. The differences between our results and the model by Debayle & Kennett (2000a) can be due to (1) the independence of the constraints from the waveform inversion used in Method I and the decorrelation of their errors; (2) the treatment of the third dimension (3-D inversion in Method I but a stack of layers in Method II); (3) somewhat different distribution of earthquake sources; (4) different data selection criteria (e.g., the frequency bands used); (5) balancing of fundamental-mode and higher-mode information through the use of one (Method II) or more (Method I) data windows; and (6) the crustal models used to apply shallow corrections.
The advantage of Method II is that a formal resolution matrix can be computed, but this comes at a substantial computational expense. We have shown that the information contained therein is qualitatively similar to ad hoc checkerboard response tests. Furthermore, the formal resolution is 2-D since the full 3-D covariance matrices of data and model are not used in Method II. As a result, variations in lateral resolution with depth and radial smearing can not be investigated efficiently with Method II. In contrast, Method I allows us to evaluate the uncertainty as well as the radial smearing in all parts of the 3-D model.
We have obtained a new shear-wave speed model for the Australian upper mantle which includes azimuthal anisotropy. The inversion presented here was based on 30 per cent more data than used by Simons et al. (1999), but we confirm their conclusions regarding the aspherical variation of isotropic wave speed in the Australian upper mantle. In particular, the conclusion that the high correlation between the westward increase in wave speed and lithospheric age holds only for depths less than 200 km is confirmed.
We have determined the spatial variation of azimuthal anisotropy (transverse anisotropy with a horizontal axis of symmetry) and investigated the trade-off between heterogeneity and anisotropy. We conclude that with the exception of westernmost Australia, where data coverage is worst, the level of anisotropy that is produced by a combination of uneven data coverage and heterogeneity is at most 20 per cent of the magnitude of anisotropy required by the data. In general the uncertainty in the angle of fast polarization is about 20°; locally the directions are more tightly constrained, but there are also regions where the uncertainty is much larger, in particular in western Australia. In the top 150 km the inversions yield fairly strong anisotropy, and the orientation of the fast axis changes in a complex manner. In contrast to the wave speed heterogeneity, there is no obvious relationship between anisotropy and the large scale geological domains that we considered, suggesting that the length scale for variation is much smaller than the crustal elements identified by Shaw et al. (1995). The horizontal length scale over which anisotropy changes increases between 150 and 200 km depth.
In the northern central region we observe a predominance of N–S directions throughout the depth range ≥140 km. Within the substantial uncertainty on their direction, these N–S oriented domains are subparallel to the direction of absolute plate motion (DeMets et al. 1990). If shearing due to plate motion is the dominant process generating seismic anisotropy, then it is perhaps surprising that so many of the strong N–S anisotropic directions are confined to central Australia and start at depths as shallow as 140 km or even less, i.e., well within what we interpret as continental lithosphere. Furthermore, complex flow around stable continental interiors could yield anisotropic fabrics in relationships that are not straightforward to correlate with a single motion vector of the entire Indo-Australian plate. If fossil anisotropy is the dominant process, then one needs to explain the change in length scale and behavior between 150 and 200 km depth for the regions outside of the Central Australian domain. Most likely, an interpretation will involve a combination of fossil and present-day anisotropy, but modelling this will require further study.
We have used the approach suggested by Montagner et al. (2000) to predict SKS splitting times and directions from our 3-D anisotropic model. Much as the original inversion result for depths above 200 km, the predicted SKS splitting map does not reveal a predominance of directions parallel to plate motion. The first order difference between our result and observations of SKS and SKKS splitting is that the surface waves seem to prefer a higher level of horizontal transverse anisotropy than the body waves. This is true, in particular, for the Precambrianz shields where shear-wave birefringence yielded split times of less than 1 s whereas the surface-wave inversion predicts times that are at least twice that large. The amplitude of anisotropy relative to the velocity anomalies, to which splitting delay times are sensitive, is not well resolved, but even with correct delay times, there remains a significant mismatch between the predicted and the observed splitting directions. In addition, frequency-dependent effects or the presence of dipping axes of symmetry might be responsible. We argue that, at least conceptually, the observed shear-wave splitting could be reconciled with the inferences from surface-wave propagation with a model of anisotropy in which the correlation wavenumbers and magnitude of anisotropy increase towards the surface.
Acknowledgments
Funding for the instruments and field campaigns involved in the skippy and kimba experiments was provided by the Australian National University. We thank Eric Debayle, Don Forsyth, Hrafnkell Kárason, Brian Kennett, Sergei Lebedev and Albert Tarantola for stimulating discussions, Takashi Iidaka for providing his shear-wave splitting analyses and Mark Leonard and the IRIS Data Management Center for data access. Guust Nolet, Steve Ward, and an anonymous reviewer provided comments which improved the manuscript. RvdH and FJS thank the Institut de Physique du Globe in Paris, where parts of this paper were written. This research was supported by NSF (grant EAR-0001136) and the David and Lucile Packard Foundation, and it is IPGP contribution # 1793.
