Shear wave structure from joint analysis of seismic and seafloor compliance data

Summary

1 Introduction

The presence of fluid within the marine crust has profound geological and geophysical implications. Crustal accretion at mid-ocean ridges occurs by crystallization in magma chambers beneath the ridge axis, while heat is removed by hydrothermal circulation through the porous upper crust. Fluid carried down into the mantle by subducting slabs can influence seismogenic and volcanic processes, and small fluid inclusions in homogeneous rocks can lead to effective anisotropy.

The key to constraining the amount and location of fluid within the crust is the determination of shear wave structure. Fluid regions can be identified from their zero shear wave speed, while estimates of porosity and the pore aspect ratio can be made from a knowledge of the Poisson ratio (ν) and measurements of attenuation.

Marine seismic methods have been very successful at revealing compressional wave structure in high-resolution detail, but information concerning shear wave structure must be inferred from waves that have undergone a mode conversion in the seabed (since the source is within water and therefore generates compressional waves only). These converted phases are often difficult to identify (because of their low amplitude or the simultaneous arrival of other waves), and consequently shear wave structure is harder to determine.

A complementary method for constraining seafloor structure is to use seafloor compliance measurements. Long-period waves travelling across the surface of the ocean produce pressure fluctuations at the seafloor. These fluctuations induce displacements in the sea bed, and the transfer function between the displacement and the applied pressure (multiplied by the forcing wavenumber) is known as the seafloor compliance. Compliance is measured experimentally using data from a broad-band seismometer and pressure gauge placed on the seafloor. Seafloor compliance is principally sensitive to shear modulus, and especially to regions of low shear modulus associated with the presence of fluids (Crawford 1991), but has lower vertical resolution than seismic reflection methods.

Seismic and compliance techniques have complementary advantages: seismic methods are excellent at revealing high-resolution information concerning compressional wave structure, while compliance methods are most sensitive to shear wave structure. Using both types of data it should be possible to better characterize the velocity structure of the oceanic crust.

In this study we consider as an example the determination of on-axis shear wave structure at 9°33′N on the East Pacific Rise (EPR) using both seismic and compliance data. We restrict attention to the uppermost crust: layer 2A, which is generally interpreted to consist of extrusive rocks, and layer 2B, the underlying sheeted dyke complex of the volcanic structure. Some of the best available constraints on young crustal structure at a mid-ocean ridge come from active source seismic experiments performed at 9°–10°N on the EPR (e.g. Vera 1990; Kent 1993; Collier & Singh 1997); this has also been the location for several compliance measurements (Crawford 1999; Crawford & Webb 2002).

We begin by using the on-axis expanding spread profile (ESP) study of Vera (1990) to determine shear wave structure from the traveltimes of pS arrivals. These are waves that undergo a P to S conversion at the layer 2A/2B interface and subsequently turn within layer 2B (Christeson 1997). An independent model of shear wave structure is then obtained using seafloor compliance data from the same location. We conclude by considering possible explanations for the differences in results given by the two methods.

2 Seismic Modelling

We determine a seismic shear wave model using the data from the on-axis ESP5 of Vera (1990). The experimental data are shown in Fig. 1(a). We have applied a Butterworth filter with a low cut of 10 Hz, a high cut of 15 Hz and a 48 dB/octave rolloff. The filtering allows us to see more clearly the refracted arrivals: at around 4.5 s (reduced time) and 8 km offset, one can see pS arrivals with a slope that corresponds to an apparent velocity of 3.34 km s⁻¹.

Figure 1

(a) Seismic data from the (on-axis) ESP5 of Vera (1990), time reduced by 6 km s⁻¹. The pS arrivals lie parallel to and below the thick black line. (b) Synthetic seismograms calculated using the Vera et al. S-wave velocity model. (c) Synthetic seismograms calculated using our best-fitting S-wave model.

Open in new tabDownload slide

We model these data by assuming a laterally uniform velocity structure and calculating synthetic seismograms using a reflectivity method (Fuchs & Muller 1971). The details of the original P- and S-wave models of Vera et al. are given in Table 1. A combination of τ–p inversion and x–t forward modelling was used to determine their P-wave model; their S-wave model is calculated using a ratio of V_S/V_P= 0.43 (corresponding to ν= 0.39) in layer 2A and V_S/V_P= 0.54 (corresponding to ν= 0.29) in layer 2B. Densities are found from the empirical relationship ρ= 1.85 + 0.165V_P (Christensen & Shaw 1970, ρ is in units of g cm⁻³, V_P is in units of km s⁻¹).

1

The original seismic model determined by Vera (1990). The layer 2A/2B boundary is at a depth of 0.254 km, and the depth of the axial magma chamber, redetermined by Kent (1993), is 1.48 km.

Open in new tabDownload slide

Using this initial velocity model, the synthetic pS arrivals show an apparent velocity of 3.04 km s⁻¹, much slower than found in the data (the synthetics are shown in Fig. 1b). To produce pS arrivals that agree with the data, we keep the V_P values of Vera et al., and make adjustments to their S-wave model in layer 2B. We find that the model must be modified by increasing V_S throughout layer 2B, on average by around 0.35 km s⁻¹. Fig. 2 shows the original Vera et al. shear wave model together with our best-fitting model. We find a velocity discontinuity of 0.6 km s⁻¹ at the layer 2A/2B interface, and a V_S of 3.59 km s⁻¹at the base of layer 2B. Within layer 2B, ν varies between 0.22 and 0.25. The synthetic seismograms for this model are shown in Fig. 1(c).

Figure 2

Shear wave velocity models derived from seismic (dashed line) and compliance (solid line) data. The original model of Vera (1990) is shown with a dotted line.

Open in new tabDownload slide

We can check that the sloping arrivals seen on the seismograms are indeed pS waves and not reflected phases by replacing the sharp layer 2A/2B interface with a 300 m wide smooth transition zone. The resulting synthetics show no changes in the nature of the pS arrivals, and so we can conclude that these are not waves reflected from the layer 2A/2B boundary. It appears that the most important feature that is needed to fit the data is a high V_S throughout layer 2B.

3 Compliance Modelling

Seafloor compliance data have been obtained at 9°33′N on the rise axis (where the water depth is 2.54 km), for frequencies in the range 0.004–0.0225 Hz (corresponding to forcing wavelengths of 3–38 km). We again assume a 1-D sea bed model, and perform forward modelling of compliance (Crawford 1991), using the same philosophy as for the seismic modelling: we keep as fixed the P-wave model of Vera (1990), and adjust their S-wave model so as to give the best fit to the data.

Using the shear wave model of Vera et al. (constant V_S/V_P ratios of 0.43 in layer 2A and 0.54 in layer 2B) produces a poor fit to the data. The fit becomes even worse with the V_S model obtained in the previous section using the seismic data. To fit the compliance data we need shear velocities in both layers 2A and 2B, which are significantly lower than those required by the seismic data: for our best-fitting model, V_S lies between 1.68 and 3.14 km s⁻¹within layer 2B, corresponding to a Poisson ratio in the range 0.33 ≤ν≤ 0.44. This model is shown in Fig. 2, and its fit to the data is illustrated in Fig. 3.

Figure 3

The compliance data, together with calculated compliance curves for the best-fitting model (solid line) and for the model produced from the seismic data (dashed line). The vertical lines represent the uncertainty in the compliance measurement as defined in Crawford (1991).

Open in new tabDownload slide

The compliance data require neither a sharp increase in shear velocity in layer 2A, nor a large discontinuity across the layer 2A/2B boundary, both of which are features of the model obtained from the ESP seismic data. Nevertheless, the models obtained from the two different data sets share the same qualitative structure, showing a gradual increase of velocity with depth in layer 2B, and the presence of a melt body beneath.

4 Discussion

The two methods of obtaining a model of the on-axis shear wave structure have led to rather different results, with the compliance data producing a model with velocities significantly lower than those found from the seismic data. In this section we examine the reasons for this variation. So as to obtain a quantitative understanding of what might be causing the discrepancy between seismic and compliance measurements of shear wave speed, we focus attention on the simplest possible case of a uniform half-space. Although unrealistic as a model of the EPR, this situation has the advantage of being analytically tractable. The quantitative insights that we gain from analysis of this simple example can be used as a guide to what is happening in more complex situations.

The compliance over an isotropic half-space with P-wave speed α and S-wave speed β is given by (Sorrels & Goforth 1973)

1

It is straightforward to show that for a given α (and assuming that the Poisson ratio satisfies 0 ≤ν≤ 0.5), eq. (1) has a unique solution for β: in other words, a compliance measurement and a P-wave speed allow us to determine the S-wave speed.

We are especially concerned with situations in which our estimate of the P-wave speed is incorrect, thus leading to an erroneous solution for the S-wave speed. In particular, we consider the effects of anelasticity and anisotropy, the presence of either of which may lead to a wrong P-wave speed being used in the compliance calculation.

4.1 Anelastic effects

Anelasticity within a material causes the attenuation of waves propagating through it. Causality requires that velocities are frequency dependent in material that is anelastic. For a material with (intrinsic) quality factors Q_P and Q_S, assumed to be constant across the frequency range, the relationship between the P- and S-wave velocities at two different frequencies f_c and f_s is given by (Aki & Richards 1980)

2

3

If f_s and f_c are of the same order of magnitude then the change in wave speed with frequency is negligible. The large difference between seismic frequencies and compliance frequencies, however, means that this change can be substantial, even for modest amounts of attenuation.

Consider a situation in which we take f_s= 10 Hz as a typical seismic frequency and f_c= 0.015 Hz as a typical compliance frequency. If the wave speeds at seismic frequencies are α_s and β_s, then the wave speeds at compliance frequencies will be α_c and β_c as given by eqs (2) and (3). The measured compliance of a uniform half-space is found by substituting α=α_c and β=β_c into the right-hand side of eq. (1) to give

4

If we now take this compliance measurement and solve for the S-wave velocity, but using the seismic P-wave velocity α_s (this is analogous to the calculation of a shear wave structure using the compliance data and the Vera et al. P-wave model), we will obtain an ‘apparent’ shear velocity of β_a satisfying the equation

5

In Fig. 4 we show the ratio of the apparent shear velocity to the seismic shear velocity, β_a/β_s, as a function of the P-wave quality factor Q_P. (We have taken Q_S=Q_P/2.) The solution is shown for half-spaces with Poisson ratios of 0.2, 0.25 and 0.3 – the difference between these three cases is not great. It can be seen that the shear wave velocity obtained from compliance data will be at least 15 per cent slower than that obtained from seismic data if Q_P≲ 35, and at least 25 per cent slower if Q_P≲ 20.

Figure 4

The ratio of apparent shear wave speed to seismic shear wave speed, β_a/β_s, as a function of quality factor Q_P. Solutions are shown for half-spaces with Poisson’s ratios of ν= 0.2 (dotted line), ν= 0.25 (solid line) and ν= 0.3 (dashed line).

Open in new tabDownload slide

Experimental estimates of Q_P on-axis at 9°30′N find high attenuation (Q_P≈ 10–20) in layer 2A, but lower attenuation (Q_P > 70) within layer 2B (Wilcock 1992, 1995; Christeson 1994). These values are measurements of effective Q, and thus include contributions from both scattering and intrinsic attenuation. Only the intrinsic part is responsible for velocity dispersion. The extensive fracturing of the upper oceanic crust means that scattering processes are significant, but it is difficult to make quantitative estimates of the relative importance of scattering and intrinsic attenuation.

We can nonetheless investigate the sensitivity of our preferred seismic model to the degree of intrinsic attenuation by performing viscoelastic modelling using successively larger amounts of attenuation in layers 2A and 2B. This modelling shows that the arrival time and apparent wave speed of the pS arrivals are fairly insensitive to even quite large amounts of attenuation (Q_P= 30, Q_S= 15 in layer 2A, Q_P= 40, Q_S= 20 in layer 2B). As would be expected, the amplitude of the arrivals decreases with increasing attenuation. The waveforms only begin to show significant departures from the Vera et al. data set when very large amounts of attenuation are used (Q_P < 30).

Taken in combination with the experimental results, the modelling suggests that much of the effective attenuation in layer 2B could be intrinsic. It seems that a smaller proportion of the effective attenuation in layer 2A is intrinsic, although the actual Q values may well be lower than in layer 2B. Although it should be remembered that the calculation is for a half-space, the results of Fig. 4 certainly indicate that these values are sufficient to cause a substantial decrease in the apparent shear velocity inferred from compliance measurements.

The preceding analysis made the assumption that Q_S and Q_P were independent of frequency. It is well known that laboratory measurements of the attenuation of oceanic basalts (made at MHz) produce very similar values to those found in seismic field experiments (Tompkins & Christensen 2001). The reasons for this are not clear, but most explanations of attenuation suggest that it is caused by a range of viscous squirt mechanisms involving fluid transfer between small-scale cracks. As seismic wavelengths are already very much larger than the length-scale of these inclusions, there is no reason to suspect any new behaviour at the lower frequencies (longer wavelengths) associated with compliance.

4.2 Anisotropy

As well as attenuation, we can also consider the effect of anisotropy. Seismic experiments at 9°30′N on the EPR find that the upper 2 km of crust is anisotropic, with a fast direction of P-wave propagation aligned in a horizontal direction and approximately parallel to the rise axis (Dunn & Toomey 2001). Because we have used refracted waves to determine a velocity structure using the seismic data, we are most sensitive to these horizontal velocities. The different mechanism that is responsible for compliance means that it may well have a different pattern of sensitivity, and so using the horizontal P-wave velocity is possibly incorrect.

In the Appendix we derive the following expression for the compliance of a transversely isotropic half-space (with symmetry axis in the horizontal (1) direction, pointing along axis):

6

The quantities C₁₁, C₃₃, C₁₃ and C₄₄ are the anisotropic elastic parameters (analogous to the Lamé parameters for an isotropic situation) written in standard abbreviated form. The horizontal (quasi-) P-wave velocity is directly related to C₁₁, while the vertical P-wave velocity is directly related to C₃₃ (Mavko 1998). The velocity of the SV wave in both the horizontal and vertical direction is the same and depends on C₄₄.

The sensitivity of the right-hand side of eq. (6) to relative changes in the elastic parameters is illustrated in Fig. 5. The particular case shown is for departures about an isotropic solid with ν= 0.25, but calculations for a range of other situations show very little variation in behaviour.

Figure 5

The sensitivity of the compliance of an anisotropic half-space to changes in each of the four elastic parameters: C₁₁ (solid line), C₃₃ (dotted line), C₄₄ (dashed line) and C₁₃ (dot-dashed line). The parameters are varied from values that correspond to an isotropic material with a Poisson’s ratio of 0.25.

Open in new tabDownload slide

The compliance is most sensitive to changes in C₃₃ (corresponding to the vertical P-wave velocity). It is reasonably sensitive to changes in C₄₄, and least sensitive to changes in either C₁₁ (corresponding to the horizontal P-wave velocity) or C₁₃. This result means that in an anisotropic material the compliance and seismic measurements will depend principally on velocities in different directions.

To investigate further the consequences of this differing pattern of sensitivity we consider the specific instance for which

7

8

Here e is a parameter that measures the degree of anisotropy. Taking e= 1 corresponds to an isotropic material; if e < 1 then the P-wave speed in the horizontal direction is faster than that in the vertical direction. The compliance of the half-space is given by substituting these values into eq. (6). If we then solve for the shear wave velocity, assuming that the half-space is isotropic and taking the P-wave velocity to be given by the horizontal wave speed α, we find that the apparent shear wave speed β_a satisfies the equation

9

The solution of this equation for small departures from isotropy (e≈ 1) is given by

10

For typical values of the Poisson ratio, (3 − 4ν)/(8ν) ≈ 1 and so the apparent percentage change in shear wave velocity is approximately equal to the percentage of anisotropy. A fast horizontal wave speed (e < 1) corresponds to an apparent decrease in shear wave speed, just as is seen in the experimental results. Dunn & Toomey (2001) found that the degree of anisotropy was no more than 5 per cent, so while this effect is not vanishingly small, it probably less important than the effects of anelasticity.

5 Conclusions

Seismic methods and seafloor compliance techniques are two complementary tools that can be used to constrain seafloor structure. We have performed an analysis of both types of data recorded on-axis at 9°33′N on the East Pacific Rise to obtain a shear wave structure for the upper oceanic crust. At first sight the two methods produce inconsistent results. To successfully model the seismic pS arrivals seen in the on-axis ESP of Vera (1990), high shear wave speeds (a Poisson ratio in the range 0.22 ≤ν≤ 0.25) are needed within layer 2B. Compliance data from the same location, meanwhile, require shear wave velocities that are much lower (0.33 ≤ν≤ 0.44 within layer 2B).

This inconsistency can be resolved by considering the effects of anelasticity and anisotropy on seismic and compliance measurements. Anelasticity causes velocity dispersion, and the large difference between seismic frequencies (≈10 Hz) and compliance frequencies (≈0.015 Hz), means that the wave speeds seen by the two methods are significantly different. This difference is exaggerated by using a seismic P-wave structure to infer the S-wave structure from the compliance data. The presence of anisotropy, meanwhile, makes compliance and seismic methods principally sensitive to velocities in different directions. Both effects lead to a lower shear wave velocity being inferred from compliance data than from seismic data.

A quantitative estimate of these two effects, for the analytically tractable case of a half-space, shows that the effect of anelasticity is generally more important than that of anisotropy. Furthermore, taking values of intrinsic attenuation consistent with both the data set of Vera et al. and previous experimental measurements produces a decrease in apparent shear wave speed of the same order of magnitude as the difference between the velocity structures derived from the seismic and compliance data. Although the half-space model is clearly unrealistic as a model of the EPR, it would be expected that the same pattern of results would be seen for a more complex model.

These results illustrate the importance of properly incorporating the effects of attenuation (and to a lesser extent anisotropy) into the joint analysis of seismic and compliance data. The converse is that a simultaneous joint inversion of the two data sets that deals with anelasticity and anisotropy in a consistent way will provide stronger constraints on both the velocity and attenuation structure of the upper oceanic crust. Seismic and compliance methods have complementary sensitivities to P- and S-wave structures, while the difference in velocities seen by the two methods can be tied to the degree of attenuation. Furthermore, as velocity dispersion is caused by intrinsic attenuation, a joint inversion will allow the separation of this effect from scattering attenuation. Since even small amounts of liquid can hugely increase attenuation (Singh 2000), this provides a means of inferring the presence of fluid within the crust.

Seismic studies further south on the EPR (where spreading rates are faster, and the magma supply is therefore greater) have found high values of attenuation within both layers 2A and 2B (Collier & Singh 1998). A simultaneous joint inversion of seismic and compliance data from an extended range of sites along the EPR will provide information on along-axis variation in attenuation structure, and thus give important insight into the delivery of melt from the mantle, the process of hydrothermal circulation and the mechanism of crustal accretion.

Acknowledgments

We thank A.J. Harding and J.A. Orcutt for providing ESP5 data collected by R.V. Robert D. Conrad and R.V. Thomas Washington during 1985, and are grateful for the helpful comments of an anonymous referee. TH is supported by an Entente Cordiale scholarship from the Ambassade de France, London. IPGP contribution 1924.

References

Appendix

Appendix: the Compliance of An Anisotropic Half-Space

We restrict attention to the case of transverse isotropy, in which the stress tensor depends on five independent elastic constants: in standard abbreviated notation, and for a symmetry axis in the 1-direction, these are C₁₁, C₃₃, C₁₂, C₁₃ and C₄₄. Taking the symmetry axis to be horizontal and directed along axis corresponds to the situation for the upper crust near a mid-ocean ridge (Dunn & Toomey 2001). We therefore associate x₁≡x with the horizontal coordinate and x₃≡z with depth.

The compliance is defined as the transfer function between the normal stress applied to the half-space and the normal displacement that this induces at the free surface (Sorrels & Goforth 1973; Crawford 1991). Assuming a plane-wave dependence of e^{i (k x - ωt)} (in which k is the horizontal wavenumber), the plane-strain equation of motion can be written in the ‘propagator’ form (Woodhouse 1974; Kennett 1990)

(A1)

where the stress-displacement vector b is given by

(A2)

and the matrix A by

(A3)

Note that since this is a plane-strain problem, there is no dependence on C₁₂.

The Jordan decomposition of A,

(A4)

represents a decomposition of the wavefield into upgoing and downgoing parts. That is to say that Λ and M can be partitioned as

(A5)

where exp(Λ_↑z) → 0 as z→−∞ and exp(Λ_↓z) → 0 as z→∞. The matrix M converts the stress-displacement vector into amounts of upgoing and downgoing waves (Kennett 1983).

In our case we have downgoing waves only, since the wavefield must decay as z→∞. This implies that

(A6)

and we find an expression for the compliance by analysing what this means for the relationship between stress and displacement. We introduce an admittance matrix, Δ, defined by

(A7)

and substitution of this definition into eq. (A6) gives

(A8)

Since τ_xz vanishes at the free surface, the compliance is found by multiplying Δ_2,2 by the wavenumber k and taking the quasi-static limit ω→ 0.

The preceding calculations are algebraically involved, but after some simplification (under the assumption C₁₁C₃₃ > C²₁₃) the result that emerges is

(A9)

Substituting the isotropic values of C₁₁=C₃₃=λ+ 2μ, C₁₃=λ and C₄₄=μ back into this equation gives

(A10)

(A11)

in agreement with the isotropic result.