Search

Understanding the crustal structure and the storage and movement of fluids beneath a volcano is necessary for characterising volcanic hazard, geothermal prospects and potential mineral resources. This study uses local earthquake traveltime tomography to image the seismic velocity structure beneath Nabro, an off-rift volcano located within the central part of the Danakil microplate near the Ethiopia-Eritrea border. Nabro underwent its first historically-documented eruption in June 2011, thereby providing an opportunity to analyse its post-eruptive state by mapping subsurface fluid distributions. We use a catalogue of earthquakes detected using machine learning methods to simultaneously relocate the seismicity and invert for the three-dimensional P- and S-wave velocity structures (VP, VS) and the ratio between them (VP/VS). Overall, our model shows higher than average P- and S-wave velocities, suggesting the presence of older consolidated volcanic deposits or intrusive magmatic rocks in the crust. We identify an aseismic region of low VP, low VS, and high VP/VS ratio at depths of 6-10 km b.s.l., interpreted as the primary melt storage region that fed the 2011 eruption. Above this is a zone of high VS, low VP, and low VP/VS ratio, representing an intrusive complex of fractured rocks partially-saturated with over-pressurised gases. Our observations identify the persistence of magma in the subsurface following the eruption, and track the degassing of this melt through the crust to the surface. The presence of volatiles and high temperatures within the shallow crust indicate that Nabro is a viable candidate for geothermal exploration.

Fiber-optic distributed acoustic sensing (DAS) cables are now used to monitor microseismicity during hydraulic-fracture stimulations of unconventional gas reservoirs. Unlike geophone arrays, DAS systems are sensitive to uniaxial strain or strain rate along the fiber direction and thus provide a 1C recording, which makes identifying the directionality and polarization of incoming waves difficult. Using synthetic examples, we have shown some fundamental characteristics of microseismic recordings on DAS systems for purposes of hydraulic fracture monitoring in a horizontal well in anisotropic (vertical transverse isotropy [VTI]) shales. We determine that SH arrivals dominate the recorded signals because their polarization is aligned along the horizontal cable at the near offset, although SV will typically dominate for events directly above or below the array. The amplitude of coherent shear-wave (S-wave) arrivals along the cable exhibits a characteristic pattern with bimodal peaks, the width of which relates to the distance of the event from the cable. Furthermore, we find that S-wave splitting recorded on DAS systems can be used to infer the inclination of the incoming waves, overcoming a current limitation of event locations that have constrained events to lie in a horizontal plane. Low-amplitude SV arrivals suggest an event depth similar to that of the DAS cable. Conversely, steep arrivals produce higher amplitude SV-waves, with S-wave splitting increasing with offset along the cable. Finally, we determine how polarity reversals observed in the P and SH phases can be used to provide strong constraints on the source mechanisms.

Icequakes, microseismic earthquakes at glaciers, offer insights into the dynamics of ice sheets. For the first time in the Antarctic, we explore the use of fiber optic cables as Distributed Acoustic Sensors (DAS) as a new approach for monitoring basal icequakes. We present the use of DAS for studying icequakes as a case study for the application of DAS to microseismic datasets in other geological settings. Fiber was deployed on the ice surface at Rutford Ice Stream in two different configurations. We compare the performance of DAS with a conventional geophone network for: microseismic detection and location; resolving source and noise spectra; source mechanism inversion; and measuring anisotropic shear-wave splitting parameters. Both DAS array geometries detect fewer events than the geophone array. However, DAS is superior to geophones for recording the microseism signal, suggesting the applicability of DAS for ambient noise interferometry. We also present the first full-waveform source mechanism inversions using DAS anywhere, successfully showing the horizontal stick-slip nature of the icequakes. In addition, we develop an approach to use a 2D DAS array geometry as an effective multi-component sensor capable of accurately characterising shear-wave splitting due to the anisotropic ice fabric. Although our observations originate from a glacial environment, the methodology and implications of this work are relevant for employing DAS in other microseismic environments.

Whether or not plate tectonic processes operated on a younger, hotter Earth remains ambiguous. Seismic data from new networks in the Hudson Bay region of the Canadian Shield, where the Precambrian geological record spans more than 2 billion years, offer fresh scope to address this problem. Using receiver function analyses we show that the crust of the Rae domain, which exhibits ages of Paleo- to Neoarchean (3.9–2.7 Ga), is likely felsic-to-intermediate in composition (average Vp/Vs < 1.73) and seismically transparent with a sharp Moho. There is little evidence for modern-style plate tectonics, and based on the simplicity and spatial extent of the felsic crust, models favouring vertical tectonic processes such as crustal delamination or plume activity appear better suited to the results. Data from the Hearne domain, which exhibits widespread ~ 2.7 Ga granite-and-greenstone geology, show a more complex crust with higher Vp/Vs ratios, consistent with a greater mafic component. The Trans-Hudson Orogen (THO), proposed to be a Himalayan-scale mountain belt during the Paleoproterozoic, is thought to have formed during the ~ 1.8 Ga collision of the Superior and Churchill plates. Results from the Quebec–Baffin Island segment of the THO appear to map out the first-order shape of the underthrusting Superior plate, with elevated Vp/Vs ratios likely representing the rifted margin of the Superior craton. Consistently thicker crust is observed beneath central and southern Baffin Island (~43 km), coincident with widespread high-grade metamorphic surface geology. These features can be explained by crustal thickening due to stacking of accreted terranes during continent–continent collision, analogous to the present-day Tibetan Plateau, followed by erosion. When reviewed in light of age and compositional constraints from the geological record, our seismic observations point towards secular crustal evolution from non-plate tectonic during the Paleo- to Mesoarchean evolving towards fully-developed modern-style plate tectonics during the Paleoproterozoic.

SUMMARY In situ stress monitoring of crustal rocks is desirable as it yields insights into earthquake mechanisms, volcanic eruptions and changes in hydrocarbon reservoirs. Shear wave splitting, induced by stress-controlled cracks in the shallow crust, provides a way to infer this stress. Temporal variations in these observations can be difficult to quantify due to scatter in the data and discontinuous observations. Here we present evidence of temporal variations in shear wave splitting from a continuous time-series with a high occurrence of microseismic events recorded in a borehole. We interpret these observations in terms of stress-controlled cracks and are able to infer changes in stress and, via modelling, suggest the cause of the anisotropy. Possible origins of the temporal variation in per cent anisotropy are tidal or oil-field production processes. Our results suggest that shear wave splitting is a viable probe for inferring changes in crustal stress in cracked rock.

The delineation and characterization of fracturing is important in the successful exploitation of many hydrocarbon reservoirs. Such fracturing often occurs in preferentially aligned sets; if the fractures are of subseismic scale, this may result in seismic anisotropy. Thus, measurements of anisotropy from seismic data may be used to delineate fracture patterns and investigate their properties. Here fracture‐induced anisotropy is investigated in the Valhall field, which lies in the Norwegian sector of the North Sea. This field is a chalk reservoir with good porosity but variable permeability, where fractures may significantly impact production, e.g., during waterflooding. To investigate the nature of fracturing in this reservoir, P‐wave amplitude variation with offset and azimuth (AVOA) is analyzed in a 3D ocean‐bottom cable (OBC) data set. In general, 3D ocean‐bottom seismic (OBS) acquisition leads to patchy coverage in offset and azimuth, and this must be addressed when considering such data. To overcome this challenge and others associated with 3D OBS acquisition, a new method for processing and analysis is presented. For example, a surface fitting approach, which involves analyzing azimuthal variations in AVO gradients, is used to estimate the orientation and magnitude of the fracture‐induced anisotropy. This approach is also more widely applicable to offset‐azimuth analysis of other attributes (e.g., traveltimes) and any data set where there has been true 3D data acquisition, land or marine. Using this new methodology, we derive high‐resolution maps of P‐wave anisotropy from the AVOA analysis for the top‐chalk reflection at Valhall. These anisotropy maps show coherent but laterally varying trends. Synthetic AVOA modeling, using effective medium models, indicates that if this anisotropy is from aligned fracturing, the fractures are likely liquid filled with small aspect ratios and the fracture density must be high. Furthermore, we show that the fracture‐normal direction is parallel to the direction of most positive AVO gradient. In other situations the reverse can be true, i.e., the fracture‐normal direction can be parallel to the direction of the most negative AVO gradient. Effective medium modeling or comparisons with anisotropy estimates from other approaches (e.g., azimuthal variations in velocity) must therefore be used to resolve this ambiguity. The inferred fracture orientations and anisotropy magnitudes show a degree of correlation with the positions and alignments of larger scale faults, which are estimated from 3D coherency analysis. Overall, this work demonstrates that significant insight may be gained into the alignment and character of fracturing and the stress field variations throughout a field using this high‐resolution AVOA method.

When synthetic seismograms are computed for complicated velocity models using asymptotic methods based on ray theory, the caustic and/or pseudo-caustic geometrical catastrophes that often occur in the ray manifolds can cause many difficulties, because ray amplitudes can become infinite, inaccurate or vanish in such situations. The longer the propagation path relative to the scale of the velocity inhomogeneities and the more the propagation tends to be oblique rather than normal to velocity variations, the more troublesome the catastrophes become. Although the Maslov integral method (MI) largely cures the problems of the classical geometric ray theory (GRT) at caustics (where nearby rays in geometrical space cross) and in geometric shadow zones, it suffers from serious artefacts at pseudo-caustics (the crossing points of rays drawn in the mixed geometrical and slowness space of the Maslov integral). Some methods have been devised to overcome this limitation, but none is satisfactory when it is necessary to compute automatically large suites of seismograms for complex models without continuous intervention from a skilled seismologist. In this paper, we demonstrate how combining the Maslov integral method with a stage of Kirchhoff integration can greatly reduce MI's pseudocaustic problems while retaining MI's ability to overcome the caustic and shadow problems of GRT. In many cases, the additional computational burden increases by one order of magnitude but is still much less than that required in a finite-difference modelling, and is not great for a common workstation. The Maslov–Kirchhoff (MK) technique is shown to be substationally more accurate than ordinary Maslov integration when pseudo-caustics are present. Furthermore, it is more robust and automated than the phase-partitioning technique in curing the problem of joint caustics and pseudo-caustics. We also show that the MK scheme can be more efficient than other Kirchhoff modelling techniques where the kernel of integration is obtained by GRT.

SUMMARY The Maslov asymptotic method addresses some of the problems with standard ray theory, such as caustics and shadows. However, it has been applied relatively little, perhaps because its accuracy remains untested. In this study we examine Maslov integral (MI) seismograms by comparing them with finite-diVerence (FD) seismograms for several cases of interest, such as (1) velocity gradients generating traveltime triplications and shadows, (2) wave-front bending, kinking and folding in a low-velocity waveguide, and (3) wavefield propagation perturbed by a high-velocity slab. The results show that many features of high- and intermediate-frequency waveforms are reliably predicted by Maslov’s technique, but also that it is far less reliable and even fails for low frequencies. The terms ‘high’ and ‘low’ are model-dependent, but we mean the range over which it is sensible to discuss signals associated with identifiable wave fronts and local (if complicated) eVects that potentially can be unravelled in interpretation. Of the high- and intermediate-frequency wave components, those wavefront anomalies due to wave-front bending, kinking, folding or rapid ray divergence can be accurately given by MI. True diVractions due to secondary wave-front sections are theoretically not included in Maslov theory (as they require true diVracted rays), but in practice they can often be satisfactorily predicted. This occurs roughly within a wavelength of the truncated geometrical wave front, where such diVractions are most important since their amplitudes may still be as large as half that on the geometrical wave front itself. Outside this region MI is inaccurate (although then the diVractions are usually small). Thus waveforms of high and intermediate frequencies are essentially controlled by classical wave-front geometry. Our results also show that the accuracy of MI can be improved by rotating the Maslov integration axis so that the nearest wave-front anomaly is adequately ‘sampled’. This rotation can be performed after ray tracing and it can serve to avoid pseudocaustics by using it in conjunction with the phase-partitioning approach. The eVort needed in phase partitioning has been reduced by using an interactive graphics technique. It is diYcult to formulate a general rule prescribing the limitations of MI accuracy because of model dependency. However, our experiences indicate that two space- and two timescales need to be considered. These are the pulse width in space, the length scale over which the instantaneous wave-front curvature changes, and the timescales of pulse width and significant features in the ray traveltime curve. It seems, from our examples, that when these scales are comparable, the Maslov method gives very acceptable results.

In order to test the hypothesis that seismic anisotropy in the lowermost mantle is caused by the development of a post-perovskite lattice preferred orientation, and that anisotropy can thus be used as a probe of the dynamics of the mantle's lower boundary layer, an integrated model of texture generation in D$^{\prime\prime}$ is developed. This is used to predict the elastic anisotropy of the lowermost mantle as probed by global anisotropic tomographic inversions. The model combines the current 3D mantle flow field with simulations of the deformation of post-perovskite polycrystalline aggregates. Different descriptions of single crystal plasticity can lead to model results which are anti-correlated to each other. In models where post-perovskite deformation is accommodated by dislocations moving on (010) or (100), patterns of anisotropy are approximately correlated with the results of tomographic inversions. On the other hand, in models where dislocations move on (001) patterns of anisotropy are nearly anti-correlated with tomographic inversions. If all the seismic anisotropy extracted from global anisotropic inversions is due to the presence of a lattice preferred orientation in post-perovskite in the lowermost mantle, and if the results of the tomographic inversions are not strongly biased by the sampling geometries, these results suggest that, in contrast to ideas based on the 1D anisotropic signal, deformation of post-perovskite in the lowermost mantle may be accommodated by dislocations moving on (010) or (100). Alternatively, a significant portion of the anisotropic signal may be caused by mechanisms other than the alignment of post-perovskite crystals.

SUMMARY Maslov ray summation is ‘less local’ than ordinary ray theory, because the receiver waveform depends on non-Fermat or neighbouring rays and more information about the wavefront than just local Gaussian curvature. In this way, the Maslov solution is able to remain valid at caustics, where geometrical rays and corresponding stationary points of the Maslov phase coalesce. the wavefront information is expressed via the Legendre transformation, whereby the physical wavefront is represented as the envelope of a family of tangent ‘planes’ (Snell fronts). the actual form of the Snell fronts (true planes, sections of curves or surfaces, etc.) depends on the spatial coordinates used. Given a selection for the Snell fronts and Maslov phase, one can substitute the Maslov integral solution directly into the wave equation and obtain a transport equation for the Maslov amplitude. This direct substitution is analogous to that used in ordinary ray theory and avoids pseudo-differential operators. Sometimes the relative curvature of the physical wavefront and a tangential Snell front is zero. the envelope-forming process breaks down, because the local correspondence between the physical front and the Snell fronts is not one to one and invertible. This situation corresponds to a so-called ‘pseudo-caustic’ (slowness-domain caustic or telescopic point) in the Maslov solution. Pseudo-caustics are not real. A particular ray from the source may touch a pseudo-caustic at some time in one coordinate system, but in another system this ray will not have a pseudo-caustic (at the same time and place). It is easy to design a change of coordinates (e.g. from cartesian to curvilinear or polar) to deform a single-valued traveltime function appropriately, but a multi-valued or folded wavefront, as at a physical or real caustic, is less simple. Catastrophe theory is concerned with putting multi-valued functions into ‘normal forms’ which do not have psuedo-caustics. the manifold here is ‘Lagrangian’ and V. I. Arnold showed that a special type of deformation or ‘canonical transformation’ must be used. A ‘Lagrangian equivalence’ consists of a deformation of the ‘base’ (x-space) and/or the addition of a function on the base. the latter simply means factoring out an appropriate reference phase before Legendre transformation and we have found that this simple step is often sufficient for removing pseudo-caustics. It requires no new numerical work, only an inspection or understanding of the ray-tracing results at hand. We present some body-wave computations using the reference-phase technique for models with real caustics in 2-D and for a single-valued wavefront in 3-D. We point out that a Lagrangian equivalence may be used to turn a maximum of the Maslov phase function into a minimum. This has no effect on the frequency-domain solution, but may affect the causality of the computed waveform when the Chapman method is used to obtain the time-domain response. Causality is a property which one may need to impose explicitly. Only the non-delta or one-sided function part of the response (waveform tail) is affected by this consideration. Although zeroth-order Maslov theory correctly describes the severe waveform is clear from Secdistortion due to wavefront catastrophes, it may not adequately model the more subtle effects of smooth wavefront bending. Zeroth-order Maslov theory contains some but not all of the first-order (ω−1) terms of ordinary asymptotic ray theory. First-order Maslov theory is needed for complete consistency up to ω−1. Experimentation will several different zeroth-order Maslov representations is a simple, rapid way to ascertain the potential importance of thse more subtle waveform effects. If the waveform tails are too strong, the assumption that the Maslov (and ray theory) amplitude function can be expanded in powers of ω− may break down. Numerical integration of a wave equation is then necessary.

Summary Traveltimes and amplitudes for P-waves from a 500 km-deep source in the Kuril subduction zone have been synthesized by ray tracing in smooth 3-D models that allow general anisotropy and inhomogeneity. the aim is to compare the effects of proposed anisotropy in or near slabs with those of lateral heterogeneity alone. to concentrate on these effects, the source position, slab thickness (90 km), dip (63°) and velocity anomaly (5 per cent) are held constant. Results are presented for isotropic models with slab penetration to 670km and 1000 km. Anisotropic models with 670 km-deep slabs have anisotropy within the slab (Anderson 1987) and in a 10° wedge above the slab (Ribe 1989a; McKenzie 1979). the resulting wavefront topology is never as simple as that in a laterally homogeneous reference Earth and there is strong model dependence of shadow zones, caustics and areas of multipathing. Rays are traced through slab models defined by 3-D cubic-spline interpolation of up to 21 elastic constants. Outside the slab region, 1-D ray tracing through PREM and spherical trigonometry are used to complete the ray path. the results illustrate the importance of using both traveltime and amplitude information when interpreting slab structure from teleseismic data. Some anisotropic slab models have been found which produce large (>2 s) traveltime residuals that are similar in many parts of the world to those for the deep isotropic model, but the amplitude patterns are substantially different. the model with a deep isotropic slab produces a narrow band of large traveltime residuals (3 s) and high amplitudes in a region across northern Canada. This feature is due to the focusing of rays that have travelled down the high-velocity core of the deep slab. Regions where ray theory fails (i.e. caustics) are obvious through multipathing and amplitude singularities. Hilbert-transforms and Airy-type decay caustics should be observed in many places if the models presented are good representations of the Earth. Multipathing in along-strike regions is a pervasive feature of the models considered and the degree of such multipathing is highly dependent on the nature of the slab-boundary velocity gradients. the model with anisotropy above the slab produces multipathing (traveltime triplication) in the down-dip region (i.e. a narrow region through Europe). Identifying such non-linear or ‘catastrophic’ features in teleseismic data is potentially more diagnostic than linearized interpretations (automated inversion). Overall, the results show that a range of conservative models representing a range of structural theories can encompass a wide range of wavefront consequences. Drs M. Weber and V. Cervený are thanked for helpful reviews of the manuscript. Advice and reprints from Dr K. Fujita are appreciated and Dr D. L. Anderson is thanked for suggesting the anisotropic slab model with no isotropic velocity anomaly. the authors acknowledge the support of NSERC (Canada) through Operating Grant number A1465. J-M.K. is supported by an Amoco Postgraduate Scholarship and an Ontario Graduate Scholarship.

SUMMARY Lighthill and others have expressed the ray-theory limit of Green’s function for a point source in a homogeneous anisotropic medium in terms of the slowness-surface Gaussian curvature. Using this form we are able to match with ray theory for inhomogeneous media so that the final solution does not depend on arbitrarily chosen ‘ray coordinates’ or ‘ray parameters’ (e.g. take-off angles at the source). The reciprocity property is clearly displayed by this ‘ray-coordinate-free’ solution. The matching can be performed straightforwardly using global Cartesian coordinates. However, the ‘ray-centred’ coordinate system (not to be confused with ‘ray coordinates’) is useful in analysing diffraction problems because it involves 2 X 2 matrices not 3 x 3 matrices. We explore ray-centred coordinates in anisotropic media and show how the usual six characteristic equations for three dimensions can be reduced to four, which in turn can be derived from a new Hamiltonian. The corresponding form of the ray-theory Green’s function is obtained. This form is applied in a companion paper.

SUMMARY An S-wavefront from an isotropic region is expected to separate into two fronts when it passes into a gradually more anisotropic region. Standard ray expansions may be used to continue the waves in the anisotropic region when these two S-wavefronts have separated sufficiently. However, just inside the anisotropic region the two S-waves interfere with an effect that is stronger than the usual ω-1 corrections of the ray method. A waveform distortion can occur and this should be considered when modelling S-waves in, e.g., subduction zones with regions of isotropy grading into regions of anisotropy. The interference is studied here by local analysis of an integral equation obtained by the Green's function method. It is found that if the elasticity and its first two derivatives are continuous at the isotropy/anisotropy border, then zeroth-order ray theory may still be used to continue the incident wave into the anisotropic region. The incident displacement is simply resolved into two definite directions at the point where the anisotropy begins. These two directions are the limits of the unique eigenvectors on the anisotropic rays as the point of isotropy (onset of splitting) is approached. If the nth derivative of the elasticity is discontinuous at the isotropy/anisotropy border, then the scattering integral which describes the interference makes a correction to ray theory which is O(ω-1/n+1) in magnitude. Hence, the interference effect is stronger when the emergence of anisotropy is more gradual. Although the corrections are given by simple expressions, it is not reasonable to specify numerical velocity models up to such high-order derivatives. For a smooth interpolation scheme, such as cubic splines, it is more practical to monitor the splitting rays obtained by ray tracing and to use the best-fitting ‘equivalent’ high-order discontinuity. This will lead to an estimate of the importance of the correction terms. An example is given for a subduction zone model involving olivine alignment in the mantle-wedge above the slab.

[1] The Seychelles plateau is a prime example of a microcontinent, yet mechanisms for its creation and evolution are poorly understood. Recently acquired teleseismic data from a deployment of 26 stations on 18 islands in the Seychelles are analyzed to study upper mantle seismic anisotropy using SKS splitting results. Strong microseismic noise is attenuated using a polarization filter. Results show significant variation in time delays (δt = 0.4–2.4 s) and smooth variations in orientation (ϕ = 15°–69°, where ϕ is the polarization of the fast shear wave). The splitting results cannot be explained by simple asthenospheric flow associated with absolute plate motions. Recent work has suggested that anisotropy measurements for oceanic islands surrounding Africa can be explained by mantle flow due to plate motion in combination with density-driven flow associated with the African superswell. Such a mechanism explains our results only if there are lateral variations in the viscosity of the mantle. It has been suggested that the Seychelles are underlain by a mantle plume. Predictions of flow-induced anisotropy from plume-lithosphere interaction can explain our results with a plume possibly impinging beneath the plateau. Finally, we consider lithospheric anisotropy associated with rifting processes that formed the Seychelles. The large variation in the magnitude of shear wave splitting over short distances suggests a shallow source of anisotropy. Fast directions align parallel to an area of transform faulting in the Amirantes. Farther from this area the orientation of anisotropy aligns in similar directions as plate motions. This supports suggestions of transpressive deformation during the opening of the Mascarene basin.

Geologists can usually interpret the rocks they encounter on Earth in the light of tectonic and volcanic processes presently operating at the plate boundaries. This approach works well for relatively young rocks (Phanerozoic: younger than 550 million years old), but for the older rocks that formed during Precambrian times (more than 550 million years old), the “plate tectonic” assumption must ultimately break down. Processes operating on the younger, hotter Earth would have been quite different to those we see today. Gathering detailed evidence preserved deep within the plates in the ancient cores of the continents (“shields”), is thus essential to our understanding of the early Earth. This can be achieved using data from dense seismograph networks, but building and maintaining them in remote areas is both logistically and financially challenging; innovative station and equipment designs are required to deliver the success enjoyed in gentler climes. The Hudson Bay Lithospheric Experiment (HuBLE), a recent UK-Canadian venture in Arctic Canada, has addressed these issues in order to place fundamental constraints on Earth structure beneath the Canadian Shield. The resulting data provide a tantalizing hint as to the processes that operated on the youthful Earth.

This paper considers only the simple part of Fourier Integral theory required for the work in hand. A special type of Henrici rolling sphere harmonic analyzer has been built which makes possible a large reduction in the amount of labor involved in dealing graphically with Fourier Integral operations. Examples of the use of this procedure as applied to several geophysical instrument problems are given.