The imprint of crustal density heterogeneities on regional seismic wave propagation

Density heterogeneities are the source of mass transport in the Earth. However, the 3-D density structure remains poorly constrained because travel times of seismic waves are only weakly sensitive to density. Inspired by recent developments in seismic waveform tomography, we investigate whether the visibility of 3-D density heterogeneities may be improved by inverting not only travel times of specific seismic phases but complete seismograms. As a first step in this direction, we perform numerical experiments to estimate the effect of 3-D crustal density heterogeneities on regional seismic wave propagation. While a finite number of numerical experiments may not capture the full range of possible scenarios, our results still indicate that realistic crustal density variations may lead to travel-time shifts of up to ∼ 1 s and amplitude variations of several tens of percent over propagation distances of ∼ 1000 km. Both amplitude and travel-time variations increase with increasing epicentral distance and increasing medium complexity, i.e. decreasing correlation length of the heterogeneities. They are practically negligible when the correlation length of the heterogeneities is much larger than the wavelength. However, when the correlation length approaches the wavelength, density-induced waveform perturbations become prominent. Recent regional-scale full-waveform inversions that resolve structure at the scale of a wavelength already reach this regime. Our numerical experiments suggest that waveform perturbations induced by realistic crustal density variations can be observed in high-quality regional seismic data. While density-induced travel-time differences will often be small, amplitude variations exceeding ±10 % are comparable to those induced by 3-D velocity structure and attenuation. While these results certainly encourage more research on the development of 3-D density tomography, they also suggest that current full-waveform inversions that use amplitude information may be biased due to the neglect of 3-D variations in density.


2002)
. Recently, Koelemeijer et al. (2015) have shown that density estimates from previously used normal-mode data are not robust. However, with the incorporation of the latest data, these inferences can be improved significantly.
On regional scales, several authors jointly inverted body wave traveltimes and gravity data under the assumption seismic velocities and density are almost uniformly scaled to each other (e.g. Tondi et al., 2000Tondi et al., , 2009Maceira & Ammon, 2009).
While being correct for purely thermal density variations, this assumption prevents the detection of those interesting cases 5 where velocities and density are not simply scaled due to the presence of compositional heterogeneities.
With the steadily increasing quality of seismic data, new observables with sensitivity to 3D density variations are becoming sufficiently robust. ;  propose to use Rayleigh-wave ellipticity and local amplification measurements to estimate lithospheric density variations. The design of seismic observables with maximum sensitivity to density and minimum trade-offs to other parameters, e.g. velocities, has been suggested by Bernauer et al. (2014). 10 In addition to improving data quality, new opportunities may arise from the development of full-waveform inversion techniques that are capable of exploiting complete seismograms without being restricted to the well-known seismic phases (e.g. Chen et al., 2007;Fichtner et al., 2009;Tape et al., 2010;. As shown by Rickers et al. (2012Rickers et al. ( , 2013, the exploitation of scattered waves in full-waveform inversion can lead to substantial improvements in regional 3D velocity images. However, the potential of full-waveform inversion to better constrain density variations in the crust and upper mantle remains largely 15 unexplored.

Outline
As a first step towards full-waveform inversion for regional density structure, we present a study on the imprint of 3D density heterogeneities in the crust on seismic wave propagation in the period range from 8 -50 s. For this, we conduct a series of numerical experiments, where we analyse seismic wave propagation through random Earth models with variable complexity. 20 These models are designed to represent are range of plausible 3D heterogeneous crustal environments. While wave propagation through random media has been widely used to quantify the effect of velocity heterogeneities (e.g. Frankel & Clayton, 1986;Frankel, 1989;Igel & Gudmundsson, 1997;Furumura & Kennett, 2005;Kennett & Furumura, 2008;Meschede & Romanowicz, 2015), variations in density have so far not been considered. Our experiments are intended to (i) provide rough estimates of the amplitude and traveltime variations related to realistic density variations, and (ii) better understand the physics behind 25 density-induced waveform perturbations.
Following a presentation of the numerical setup, we will present detailed analyses of traveltime and amplitude variations induced by 3D crustal density heterogeneities. We expect scattering to be the dominant mechanism by which density heterogeneities influence the seismic signal. Scattering is most effective when scatterers are of similar size or smaller than the wavelength, which is why we will study the influence of frequency, propagation distance and medium complexity. Being fo-30 cused on a future full-waveform inversion for density, we do not consider specific seismic phases, but try to provide ensemble estimates of waveform perturbations. Given the complexity of regional-scale seismic waveforms at periods below ∼ 20 s, it is clear that this analysis can never be complete and exhaustive. It will, however, provide a first crude estimate of the impact of crustal density structure on seismic wave propagation.
2 Setup of the numerical experiments

Numerical wave propagation
To assess the impact of 3D density heterogeneities in the crust on seismic wave propagation, we compute numerical solutions 5 to the elastic wave equation (e.g. Kennett, 2001;Aki & Richards, 2002) which relates mass density ρ, the elastic tensor c ijkl , and an external force f i to the displacement field u i . With our focus being on regional wave propagation at periods below 50 s, we can safely ignore the Earth's rotation and self-gravitation. Furthermore, we restrict ourselves to an isotropic rheology. 10 For the numerical solution of equation (1), we employ the spectral-element solver SES3D (Fichtner et al., 2009;Gokhberg & Fichtner, 2016). The spectral-element method, widely used in seismological research, allows us to compute accurate numerical solutions in the presence of strong 3D heterogeneities, and without requiring special treatment of the free surface (e.g. Faccioli et al., 1996;Komatitsch & Vilotte, 1998;Komatitsch & Tromp, 1999;Peter et al., 2011).
Our computational domain is a spherical section that is ∼ 2000 by ∼ 1000 km wide and 500 km deep. As background model we 15 use the radially symmetric Preliminary Reference Earth Model (Dziewoński & Anderson, 1981), where we replace the original 24 km thick crust by a 40 km thick crust that better represents continental structure.
Since the number of receivers has no significant impact on the computational costs of the numerical simulations, we use a dense grid of 930 receivers, distributed evenly across the surface of the computational domain. In the wavefield simulations, we calculate 700 s long velocity seismograms from a strike-slip source. The complete setup is shown in figure 1.

Random media generation
With the true 3D density structure of the crust being insufficiently constrained, we use synthetic random density models in our numerical wave propagation experiments. For this, we superimpose random velocity and density variations with pre-defined correlation lengths in the horizontal and vertical directions onto the crustal part of the background model, i.e. the upper 40 km.
The spatial variations in velocity and density are statistically uncorrelated. In Appendix A, we summarise the computation of 25 3D random models based on the widely-used Fourier method.
To ensure that the amplitudes of velocity and density variations are realistic, we combine information from tomographic models and empirical velocity-density scalings. For this, we compute the root-mean square (rms) of the S velocity variations in the of these variations, and the range of different velocity-density scalings proposed in the literature. We discuss these issues in more detail in section 4.2. We note that random velocity and density models used in our simulations are on purpose spatially uncorrelated. Empirical velocity-density scalings are used only to determine plausible rms variations.

Quantification of waveform differences
In our numerical experiments, we use media with homogeneous crustal density and random 3D variations in P and S velocities 10 as reference. We then compare synthetic seismograms from the reference medium with synthetic seismograms from a medium where random variations in density are added.
Since our ultimate goal is to use complete three-component seismograms to constrain density in the Earth, we do not compare isolated and well-defined seismic phases. Instead, we compute time-and frequency-dependent traveltime and amplitude differences. For this, we bandpass filter the seismograms into a pre-defined frequency band and apply a zero-centred mov- 15 ing window w(t) that transforms a component of a seismogram u(t) into its windowed versionû τ (t) = w(t − τ )u(t). The traveltime difference δT as a function of τ is then defined as the argument of the maximum of the cross-correlation: whereû ref τ denotes the windowed seismogram for the reference medium with homogeneous crustal density. In the case of δT < 0, the wave for 3D heterogeneous density arrives earlier than the reference wave, and vice versa. Similarly, we measure 20 relative amplitude variations δA as a function of time: In the following sections, we consider three frequency bands of variable width: 0.02 -0.125 Hz (8 -50 s), 0.02 -0.067 Hz (15 -50 s), and 0.02 -0.04 Hz (25 -50 s). To stabilise the measurements, we exclude those parts of the synthetic seismograms where the average amplitude within a time window is below 5 % of the maximum within the complete trace. 25 While more information-rich quantifications of seismic waveform differences may be constructed, for instance on the basis of wavelet transforms (e.g. Kristekova et al., 2006Kristekova et al., , 2009, we prefer the traveltime and amplitude differences defined in equations (2) and (3) for their robustness and ease of interpretation. Similar quantifiers of waveform differences are frequently used in full-waveform inversion (e.g. Fichtner et al., 2008;van Leeuwen & Mulder, 2010;Bozdag et al., 2011;Rickers et al., 2012Rickers et al., , 2013. In the following sections, we present a phenomenological study on the impact of crustal density heterogeneities for media with different horizontal and vertical correlation lengths, summarised in table 2. For each experiment, we compare three-component seismograms in three different frequency bands: seismograms computed for a medium with 3D random density variations are compared with seismograms for the reference medium with homogeneous crustal density. We specifically analyse the effects 5 of frequency (section 3.2), epicentral distance (section 3.3), and medium complexity (section 3.4).

A single-receiver example
We start with the analysis of media with 200 km lateral and 20 km vertical correlation length. This will serve as a baseline for later simulations with models that have either more or less complexity. Figure 3 shows a comparison of three-component seismograms for homogeneous and heterogeneous crustal densities in the broadest frequency band from 0.020 -0.125 Hz (8 -10 50 s). Before attempting a more comprehensive analysis in the following sections, we consider a single receiver located at 910 km epicentral distance, marked by the red triangle in figure 1.
Waveform differences mostly tend to increase with increasing traveltime, in accord with the expectation that (multiply) scattered waves should arrive later than the primary waves by which they have been excited. The magnitude of the time shifts are approximately independent of the component, reaching around 0.5 s. Relative amplitude differences are largest on the E-W and 15 vertical components, where the displacement velocity itself is smallest so that the influence of scattered waves is largest. They regularly exceed 50 % in both directions, meaning that amplitudes for the heterogeneous density crust can be both twice and half as large as for the medium with homogeneous crustal density. On the N-S component, where the displacement velocity is largest, relative amplitude differences vary between ±10 %.  In addition to the dependence on the random velocity and density structure, figure 4 also reveals that the waveform differences in the lower frequency band from 0.02 -0.04 Hz are on average smaller than for the higher frequency band from 0.020 -0.125 Hz. In the following section, we will investigate this frequency dependence in more detail. This preliminary, and mostly visual, analysis of frequency dependence indicates that waveform differences are primarily caused by scattering that transfers energy from the large-amplitude N-S component onto the smaller-amplitude E-W and vertical components. Constructive and destructive interference between primary and scattered waves may cause the wave amplitudes to deviate in both directions. An increase of amplitudes may be further supported by additional wave focusing induced by 3D density heterogeneities. The approximate frequency independence of traveltime differences, however, can hardly be explained 5 with basic wave propagation intuition.

The effect of frequency
To make our analysis more comprehensive and efficient, we compute histograms of time shifts and relative amplitude differences for all 930 stations in the receiver grid. In line with our future goal, which is to use full-waveform inversion to constrain 3D density variations, we do not consider specific seismic phases, but longer time series that comprise body, surface and scattered waves. For this, we perform measurements in a 300 s long time interval, starting at the P wave arrival. This procedure To investigate this phenomenon further, we show histograms for the lowest and highest frequency bands for a single random medium realisation and for the three different components in the top row of figure 7. We observe a distinct tail of reduced time shifts of ∼ 1 s on the E-W component and for the lowest frequency band. A similar observation can be made for two 25 out of five random media, suggesting that these time shifts are not highly unlikely artefacts of an unusual random medium realisation. The bottom row of figure 7 shows a pair of synthetic seismograms at an epicentral distance of ∼ 1300 km that contributes to this tail of negative time shifts. In the lower frequency band from 0.02 -0.04 Hz waves for the 3D heterogeneous crust arrive early by ∼ 1 s. Histograms for larger epicentral distance have much bigger spread, and the isolated tail of negatve timeshifts appears only for far-away stations. We examined this distance dependence further, using synthetic data from one of epicentral distance between 100 and 300 km. A visual waveform comparison suggests that the traveltime differences may be a finite-frequency effect, meaning that waveform (amplitude) differences within short time intervals translate into time shifts when these are measured by cross-correlation within a finite frequency band. We discuss this aspect in more detail in section 4.1.

The effect of epicentral distance 5
To investigate whether density-related traveltime and amplitude differences are only local effects or accumulate with propagation distance, we plot histograms for stations in two different epicentral distance ranges: 100 -300 km and 1000 -1200 km.
We This indicates that waveform differences due to crustal density heterogeneities indeed accumulate with increasing epicentral distance, which is an essential prerequisite for the use of tomographic methods.

The effect of medium complexity
In order to reveal the physical origin of the waveform differences, we consider random media with different lateral correlation The histograms in figure 9 indicate that waveform differences induced by 3D density variations occur mostly due to scattering which is most effective when heterogeneities are equal or smaller in size than the wavelength. In a medium with 50 km lateral correlation length, the size of heterogeneities is comparable to the wavelength of waves with a maximum frequency of 0.125 Hz, and scattering becomes the dominant mechanism to perturb the wavefield. In the smooth medium with 1000 km lateral 25 correlation length, the dominant mechanism is transmission, which clearly has no significant impact on either traveltimes or amplitudes.

Seismic signatures of crustal density heterogeneities
Our numerical experiments show that 3D crustal density heterogeneities may lead to both positive and negative variations in the traveltimes and amplitudes of seismic waves. This indicates that 3D density structure leaves an imprint on regional seismic wavefields that goes beyond simple scattering attenuation of the main arrivals. 5 To understand effects which play a major role in wave propagation, we look first at the misfit histograms for different frequency bands (section 3.2). Intuitively, the histogram for the frequency band in which more misfits are accumulated should have smaller zero peak and bigger spread. However, as can be seen in figure 6, while the peak around zero is smaller for the lower frequency band for both time shifts and amplitude differences, the spread is comparable for time shifts and larger for the higher frequency band for amplitude differences. This may mean that we observe two separate effects affecting wave propagation that have 10 different impact on the observed seismograms for different frequency bands. One of those noticeable effects could be that, with increasing frequency, scattering off heterogeneities becomes less significant in favour of transmission effects. That essentially means approaching the range of the infinite-bandwidth approximation of ray theory validity, in which we lose sensitivity to density completely. Smaller scattering in higher frequencies may cause the histograms for the higher frequency band to have larger zero peaks. Contrarily, an increase in frequency corresponds to a proportional increase in relative propagation distance, 15 which means that we will observe misfits accumulated for bigger number of wavelenghts. The larger propagation distance could be then the reason behind the broader histogram spread for the higher frequency band. This effect is equivalent to moving further from the source, which is consistent with section 3.3, where we show that density-related misfits accumulate with distance. The increase of waveform differences with increasing epicentral distance also suggests that the impact of density structure is not merely a local effect, but rather an integral over the complete wavepath -an essential prerequisite for performing 20 tomography. While the results for different frequency bands are governed both by the scattering to transmission ratio and the propagation distance effect, and the results for different epicentral distances by the propagation distance effect only, for various media complexities we observe how important scattering is for sensitivity to density. The noticeable change in histogram shape in figure 9 is caused by much bigger amount of scattering for more complex media. Especially the nearly complete absence of waveform differences for long-wavelength density heterogeneities indicates that the scattering is the dominant mechanism to 25 produce these waveform differences. The energy transfer between the different components and towards later arrivals is also scattering-related.
As we show, the behaviour of density-induced misfits can be most often explained by an interplay of two effects: scattering and propagation distance. However, not all of the observed features can be interpreted on this ground. For instance, traveltime variations do not seem to exhibit a pronounced frequency-dependence, in contrast to amplitude variations that decay rapidly finite frequency band are controlled by the complex interference of direct and scattered waves, which may lead to seemingly paradoxical effects (e.g. Tong et al., 1998;Marquering et al., 1999;Dahlen et al., 2000). This may include the large negative traveltime differences shown in figure 7. Figure 10 illustrates finite-frequency traveltime changes for a 2D P-SV wavefield simulation. While interacting with a density anomaly, the wavefield undergoes reflections and amplitude changes. The actual onset time of the wavefront remains largely 5 unaffected. However, the waveform differences translate into a cross-correlation traveltime shift of ∼ 1 s.
The density-induced waveform differences that we found in our numerical experiments are above the noise level of many of today's regional-scale seismic recordings. While this indicates that density heterogeneities do leave a measurable imprint, it does not automatically imply that crustal density structure can be easily recovered in a tomographic inversion. Trade-offs with P and S velocity structure, for instance, may prevent the unambiguous reconstruction of density heterogeneities. The resolvability 10 of density structure may be analysed using principal-component analysis of finite-frequency kernels (Sieminski et al., 2009), and it may be improved by the construction of targeted misfit functionals (Backus & Gilbert, 1968Bernauer et al., 2014).

Random models of plausible Earth structure
In the absence of detailed information on crustal density structure on regional scales, we base our numerical experiments on For simplicity, we assume that the amplitude spectrum of the crustal velocity variations is white, meaning that velocity and density variations have nearly identical power at all scales considered in this study, i.e. from 50 -1000 km. The resulting velocity and density variations may be too large or too small by several percent, depending on whether a specific region is 30 long-term stable and subject to recent tectonic activity.
Uncertainties in velocity-density scalings are mostly caused by the natural scatter of the velocity-density relation in natural rocks. While S velocity typically varies less than ∼ 10 % for a given crustal P velocity, density can easily vary by more than 50 % for a given v P /v S ratio (Brocher, 2005). Despite the natural scatter, published velocity-density relationships for the continental crust as a whole show good agreement with their respective range of validity (e.g. Ludwig et al., 1970;Gardner et al., 1974;Christensen & Mooney, 1995), suggesting that the choice of a particular one does not introduce a significant bias.
In the light of these uncertainties, it must be kept in mind that the waveform variations resulting from our synthetic random 5 density heterogeneities represent a first rough estimate. It is intended to reveal the first-order effects but not the smaller details that certainly depend on the characteristics of a specific region.

Velocity bias estimation
The shifts in travel time observed here as a result of density structure may cause a bias in velocity structure obtained in tomographic models. In order to obtain an estimate of these velocity biases, we take a simplified approach. We consider the 10 highest frequency band from 0.02 -0.125 Hz and the reference medium with 200 km lateral and 20 km vertical correlation length, previously used in section 3.3 on the effect of epicentral distance. Since we do not analyse specific seismic phases, we take the time shift variances for different epicentral distances as representative values of time shifts induced by 3D density structure. Furthermore, we assume that sensitivity to velocity structure is concentrated on the great circle connecting source and receiver. 15 Based on these simplifications, we estimate that the shear velocity bias for an epicentral distance of 200 km is ∼ 25 m/s or ∼ 0.78 % relative to the upper-crustal shear velocity of PREM (Dziewoński & Anderson, 1981). For epicentral distances of 1000 km, the same bias is ∼ 21 m/s or ∼ 0.66 %. It follows that the velocity biases induced by the neglect of 3D density variations are small compared to the crustal velocity variations inferred from traveltime tomography, which are on the order of 10 %. However, depending on the tomographic resolution and data quality, the biases may be larger than the error bars.

Attenuation bias estimation
To quantify potential biases in attenuation induced by unknown 3D density structure, we adopt similar simplifications as in the previous section. In the ray theory approximation, relative amplitude differences between attenuated and attenuation-free waves are given by 25 with the background attenuation or inverse quality factor denoted by q, the epicentral distance by r, frequency by f , and velocity by v. The presence of 3D density heterogeneity induces additional relative amplitude differences, δA, that translate into apparent variations in attenuation δq,  (4) and (5), yields the fractional apparent attenuation bias δ ln q as a function of δA, δ ln q = ln(A + δA + 1) ln(A + 1) − 1 .
Equation (6) reveals that attenuation biases for a given δA have a dependence on the background attenuation q through A = e −πrf q v . The bias is larger for larger attenuation. Figure 11 shows fractional attenuation biases δ ln q as a function of the background attenuation q for an epicentral distance of 1000 km, the highest frequency band from 0.02 -0.125, and the medium Taking the variance of the relative amplitude differences of 0.11 as a representative value (see section 3.3), the apparent variations in attenuation, δ ln q, approximately range between ±20 % and ±60 % for q between 0.001 and 0.010 (Q between 100 and 1000). It follows that the apparent variations in attenuation induced by 3D density structure can be on the order of 10 tens of percent, thus being comparable to attenuation heterogeneities found on regional and global scales (e.g. Mitchell, 1995;Romanowicz, 1995;Haberland & Rietbrock, 2001;Selby & Woodhouse, 2002;Dalton et al., 2008;Kennett & Abdullah, 2011;Zhu et al., 2015). This result highlights that 3D density structure should be taken into account when using full-waveform techniques to invert for 3D variations in attenuation.

15
We presented a series of numerical experiments to study the effect of 3D crustal density heterogeneities on regional seismic wave propagation in the frequency range from 0.02 -0.125 Hz (8 -50 s). These were intended to (i) reveal the extent to which density-induced waveform perturbation may be measurable, and (ii) facilitate a better intuitive understanding of the underlying wave propagation physics. Our most important finding is that waveform perturbations induced by realistic crustal density variations can certainly be observed in modern, high-quality regional seismic data. While traveltime differences of typically less than 1 s will often be small compared to traveltime differences caused by velocity heterogeneities, amplitude variations of more than 10 % are comparable with those induced by 3D velocity structure and attenuation. This implies, on the one hand, that density structure may to some extend be constrained in future full-waveform inversions. On the other hand it suggests that current full-waveform 5 inversions that use amplitude information may be biased due to the neglect of 3D variations in density.
Appendix A: Random model generation We generate random media with the Fourier method, widely used in seismological research (e.g. Frankel & Clayton, 1986;Frankel, 1989;Igel & Gudmundsson, 1997;Klimeš, 2002;Jahnke et al., 2008), and recently extended to non-stationary and anisotropic media by Meschede & Romanowicz (2015). For this, we generate a random, uniformly distributed phase spectrum 10 ϕk between −π and π, with k being the 3D wavenumber. The random spectrum e i ϕk is then modulated by a positive, realvalued filter f (k) to yield the wavenumber-domain random model The filter f (k) is designed such that wavenumber components k i above a given threshold k i,max = 2π/λ i,min are excluded.
Computing the inverse Fourier transform, yields the space-domain random model By design, the 3D field m(x) only contains wavelengths above λ i,min in the i-direction. Finally, the random realisation m(x) is appropriately scaled and assigned to a specific medium parameter, such as density, P-or S-velocity.