Geophysical Journal International

SummaryTectonic slivers form in the overriding plate in regions of oblique subduction. The inner boundaries of the sliver are often poorly defined and can consist of well-defined faults, rotating blocks or diffuse fault systems, which pass through or near the volcanic arc. The Guanacaste Volcanic Arc Sliver (GVAS) as defined by Montero et al. (2017), is a segment of the Central American Forearc Sliver, whose inner boundary is the ∼87 km long Haciendas-Chiripa Fault System (HCFS), which is located ∼10 km behind the volcanic arc and consists of strike slip faults and pull apart steps. We characterise the current ground motion on this boundary by combining earthquake locations and focal mechanisms of the 2016 Bijagua earthquake sequence, with the surface ground deformation obtained from Interferometric Synthetic Aperture Radar (InSAR) images from the ALOS-2 satellite. The coseismic stack of interferograms show ∼6 cm of displacement towards the line of sight of the satellite between the Caño Negro fault and the Upala fault, indicating uplift or SE horizontal surface displacement. The largest recorded earthquake of the sequence was Mw 5.4, and the observed deformation is one of the smallest earthquakes yet detected by InSAR in the Central American region. Forward and inverse models show the surface deformation can be partially explained by slip on a single fault, but it can be better explained by slip along two faults linked at depth. The best-fit model consists of 0.33 m of right lateral slip on the Caño Negro fault and 0.35 m of reverse slip on the Upala fault, forming a positive flower structure. As no reverse seismicity was recorded, we infer the slip on the Upala fault occurred aseismically. Observations of the Bijagua earthquake sequence suggests the forearc sliver boundary is a complex and diffuse fault system. There are localised zones of transpression and trantenssion and areas where there is no surface expression suggesting the fault system is not yet mature. Although aseismic slip is common on subduction interfaces and mature strike-slip faults, this is the first study to document aseismic slip on a continental tectonic sliver boundary fault.

SummaryStation noise levels play a fundamental limitation in our ability to detect seismic signals. These noise levels are frequency-dependent and arise from a number of physically different drivers. At periods greater than 100 s, station noise levels are often limited by the self-noise of the instrument as well as the sensitivity of the instrument to non-seismic noise sources. Recently, station operators in the Global Seismographic Network (GSN) have deployed several Streckeisen STS-6A very broadband borehole seismometers. These sensors provide a potential replacement for the no-longer-produced Streckeisen STS-1 seismometer and the GeoTech KS-54000 borehole seismometer. Along with showing some of the initial observational improvements from installing modern very broadband seismometers at depth, we look at current limitations in the seismic resolution from earth tide periods 100,000 s (0.01 mHz) to Nyquist at most GSN sites (0.02 s or 50 Hz). Finally, we show the potential for improved observations of continuously excited horizontal Earth hum as well as the splitting of very long-period torsional modes. Both of these observations make use of the low horizontal noise levels which are obtained by installing very broadband borehole seismometers at depth.

SummaryWe develop a finite element dynamic earthquake simulator, EQsimu, to model multicycle dynamics of three-dimensional geometrically complex faults. The fault is governed by rate- and state-dependent friction (RSF). EQsimu integrates an existing finite element code EQdyna for the co-seismic dynamic rupture phase and a newly developed finite element code EQquasi for the quasi-static phases of an earthquake cycle, including nucleation, post-seismic, and inter-seismic processes. Both finite element codes are parallelized through Message Passing Interface to improve computational efficiency and capability. EQdyna and EQquasi are coupled through on-fault physical quantities of shear and normal stresses, slip-rates and state variables in RSF. The two-code scheme shows advantages in reconciling the computational challenges from different phases of an earthquake cycle, which include 1) handling time steps ranging from hundredths of a second to a fraction of a year based on a variable time stepping scheme, 2) using element size small enough to resolve the cohesive zone at rupture fronts of dynamic ruptures, and 3) solving the system of equations built up by millions of hexahedral elements.EQsimu is used to model multicycle dynamics of a 3D strike-slip fault with a bend. Complex earthquake event patterns spontaneously emerge in the simulation, and the fault demonstrates two phases in its evolution. In the first phase, there are three types of dynamic ruptures: ruptures breaking the whole fault from left to right, ruptures being halted by the bend, and ruptures breaking the whole fault from right to left. As the fault bend experiences more ruptures, the zone of stress heterogeneity near the bend widens and the earthquake sequence enters the second phase showing only repeated ruptures that break the whole fault from left to right. The two-phase behaviors of this bent fault system suggest that a 10° bend may conditionally stop dynamic ruptures at the early stage of a fault system evolution and will eventually not be able to stop ruptures as the fault system evolves. The nucleation patches are close to the velocity strengthening region. Their sizes on the two fault segments are different due to different levels of the normal stress.

SummaryTomography is one of the cornerstones of geophysics, enabling detailed spatial descriptions of otherwise invisible processes. However, due to the fundamental ill-posedness of tomography problems, the choice of parametrizations and regularizations for inversion significantly affect the result. Parametrizations for geophysical tomography typically reflect the mathematical structure of the inverse problem. We propose, instead, to parametrize the tomographic inverse problem using a geologically motivated approach. We build a model from explicit geological units that reflect the a priori knowledge of the problem. To solve the resulting large-scale nonlinear inverse problem, we employ the efficient Ensemble Kalman Inversion scheme, a highly parallelizable, iteratively regularizing optimizer that uses the ensemble Kalman filter to perform a derivative-free approximation of the general iteratively-regularized Levenberg-Marquardt method. The combination of a model specification framework that explicitly encodes geological structure and a robust, derivative-free optimizer enables the solution of complex inverse problems involving non-differentiable forward solvers and significant a priori knowledge. We illustrate the model specification framework using synthetic and real data examples of near-surface seismic tomography using the factored eikonal fast marching method as a forward solver for first arrival travel times. The geometrical and level set framework allows us to describe geophysical hypotheses in concrete terms, and then optimize and test these hypotheses, helping us to answer targeted geophysical questions.

SummaryAccurate synthetic seismic wavefields can now be computed in 3D earth models using the spectral element method (SEM), which helps improve resolution in full waveform global tomography. However, computational costs are still a challenge. These costs can be reduced by implementing a source stacking method, in which multiple earthquake sources are simultaneously triggered in only one teleseismic SEM simulation. One drawback of this approach is the perceived loss of resolution at depth, in particular because high-amplitude fundamental mode surface waves dominate the summed waveforms, without the possibility of windowing and weighting as in conventional waveform tomography. This can be addressed by redefining the cost-function and computing the cross-correlation wavefield between pairs of stations before each inversion iteration. While the Green’s function between the two stations is not reconstructed as well as in the case of ambient noise tomography, where sources are distributed more uniformly around the globe, this is not a drawback, since the same processing is applied to the 3D synthetics and to the data, and the source parameters are known to a good approximation. By doing so, we can separate time windows with large energy arrivals corresponding to fundamental mode surface waves. This opens the possibility of designing a weighting scheme to bring out the contribution of overtones and body waves. It will also make it possible to balance the contributions of frequently sampled paths versus rarely sampled ones, as in more conventional tomography. Here we present the results of proof of concept testing such an approach for a synthetic 3-component long period waveform dataset (periods longer than 60 s), computed for 273 globally distributed events in a simple toy 3D radially anisotropic upper mantle model which contains shear wave anomalies at different scales. We compare the results of inversion of 10,000 s long stacked time series, starting from a 1D model, using source stacked waveforms and station-pair cross-correlations of these stacked waveforms in the definition of the cost function. We compute the gradient and the Hessian using normal mode perturbation theory, which avoids the problem of cross-talk encountered when forming the gradient using an adjoint approach. We perform inversions with and without realistic noise added and show that the model can be recovered equally well using one or the other cost function. The proposed approach is computationally very efficient. While application to more realistic synthetic datasets is beyond the scope of this paper, as well as to real data, since that requires additional steps to account for such issues as missing data, we illustrate how this methodology can help inform first order questions such as model resolution in the presence of noise, and trade-offs between different physical parameters (anisotropy, attenuation, crustal structure etc..) that would be computationally very costly to address adequately, when using conventional full waveform tomography based on single-event wavefield computations.

SummaryStrongly asymmetric cross-correlation functions (CCF) of ambient noise indicate an inhomogeneous distribution of the noise sources, and can be used to locate these sources. With this premise, a grid-search procedure is applied to explore the frequency dependent source locations of the single- (SF) and double-frequency (DF) microseisms recorded in eastern North American margin (ENAM). The frequency dependent Rayleigh wave group velocities are determined based on CRUST1.0 model and the Preliminary Reference Earth model and then the theoretical travel time differences from each hypothetical source to all station pairs are calculated by the fast marching method. Then a misfit function at each hypothetical source is defined to take the difference between the calculated travel time differences and the time lags observed from CCFs. A factor of reliability is defined considering the Rayleigh wave attenuation, the correlation coefficients between the microseisms and the Wavewatch III hindcasts of ocean wave spectra, and the efficiency of conversion from DF energy in water to seismic Rayleigh waves as a function of water depth and frequency. Together with the correlation analyses of the power spectra densities of the microseisms among stations, the primary source regions of SF and DF microseisms recorded in ENAM are estimated: 1) the sources of SF microseisms are likely distributed in the continental shelf dominantly as well as adjacent deep ocean area in western North Atlantic Ocean, and 2) the long-period and short-period DF microseisms appear generated in deep ocean near the continental slope and in the continental shelf and slope of ENAM, respectively.

SummaryMonitoring temporal changes of volcanic interiors is important to understand magma, fluid pressurization and transport leading to eruptions. Noise–based passive seismic monitoring using coda wave interferometry is a powerful tool to detect and monitor very slight changes in the mechanical properties of volcanic edifices. However, the complexity of coda waves limits our ability to properly image localized changes in seismic properties within volcanic edifices. In this work, we apply a novel passive ballistic wave seismic monitoring approach to examine the active Piton de la Fournaise volcano (La Réunion island). Using noise correlations between two distant dense seismic arrays, we find a 2.4 % velocity increase and -0.6 % velocity decrease of Rayleigh–waves at frequency bands of 0.5–1 Hz and 1–3 Hz, respectively. We also observe a -2.2 % velocity decrease of refracted P–waves at 550 m depth at the 6–12 Hz band. We interpret the polarity differences of seismic velocity changes at different frequency bands and for different wave types as being due to strain change complexity at depth associated with subtle pressurization of the shallow magma reservoir. Our results show that velocity changes measured using ballistic waves provide complementary information to interpret temporal changes of the seismic properties within volcanic edifices.

SummaryThe method of ScS reverberation migration is based on a “common reflection point” analysis of multiple ScS reflections in the mantle transition zone (MTZ). We examine whether ray-theoretical traveltimes, slownesses, and reflection points are sufficiently accurate for estimating the thickness H of the MTZ, defined by the distance between the 410-km and 660-km phase transitions. First, we analyze ScS reverberations generated by 35 earthquakes and recorded at hundreds of seismic stations from the combined Arrays in China, Hi-NET in Japan, and the Global Seismic Network. This analysis suggests that H varies by about 30 km and therefore that dynamic processes have modified the large-scale structure of the MTZ in eastern Asia and the western Pacific region. Second, we apply the same procedure to spectral-element synthetics for PREM and two 3-D models. One 3-D model incorporates degree-20 topography on the 410 and 660 discontinuities, otherwise preserving the PREM velocity model. The other model incorporates the degree-20 velocity heterogeneity of S20RTS and leaves the 410 and 660 flat. To optimize reflection point coverage, our synthetics were computed assuming a homogeneous grid of stations using 16 events, four of which are fictional. The resolved image using PREM synthetics resembles the PREM structure and indicates that the migration approach is correct. However, ScS reverberations are not as strongly sensitive to H as predicted ray-theoretically because the migration of synthetics for a model with degree-20 topography on the 410 and 660: H varies by less than 5 km in the resolved image but 10 km in the original model. In addition, the relatively strong influence of whole-mantle shear-velocity heterogeneity is evident from the migration of synthetics for the S20RTS velocity model and the broad sensitivity kernels of ScS reverberations at a period of 15 s. A ray-theoretical approach to modeling long-period ScS traveltimes appears inaccurate, at least for continental-scale regions with relatively sparse earthquake coverage. Additional modeling and comparisons with SS precursor and receiver function results should rely on 3-D waveform simulations for a variety of structures and ultimately the implementation of full wave theory.

SummaryThe nutation harmonic terms are commonly determined from celestial pole offset series produced from VLBI time delay analysis. This approach is called an indirect approach. As VLBI observations are treated independently for every session, this approach has some deficiencies such as a lack of consistency in the geometry of the session. To tackle this problem, we propose to directly estimate nutation terms from the whole set of VLBI time delays, hereafter referred as a direct approach, in which the nutation amplitudes are taken as global parameters. This approach allows us to reduce the correlations and the formal errors and gives significant discrepancies for the amplitude of some nutation terms. This paper is also dedicated to the determination of the Earth resonance parameters, named Polar Motion, Free Core Nutation, and Free Inner Core Nutation. No statistically significant difference has been found between the estimates of resonance parameters based upon ’direct’ and ’indirect’ nutation terms. The inclusion of a complete atmospheric-oceanic non-tidal correction to the nutation amplitudes significantly affected the estimates of the Free Core Nutation and the Free Inner Core Nutation resonant frequencies. Finally, we analyzed the frequency sensitivity of Polar Motion resonance and found that this resonance is mostly determined by the prograde nutation terms of period smaller than 386 days.

SummaryKnowledge about the stress state in sedimentary basins gives insight into geodynamics of a given region, natural fracture development, and is important in design of underground engineering operations, such as hydraulic fracturing. As the direct stress measurements are expensive, usually very limited amount of data is available, and the stress state assessment bases on theoretical models. In the present work, we review the commonly used stress prediction models. We focus especially on the ones which take into account material viscoelasticity, and stress transfer between layers characterized by different mechanical properties. We extend the stress driven elastic model to material viscoelasticity, and we apply it to predict stress changes during last glacial cycle in the Baltic basin, northern Poland. We conclude that neglecting material viscoelasticity in creeping rocks like shales or rock salt may lead to erroneous stress prediction, and that coupling of the layers induces stress transfer among layers, and together with stress relaxation in ductile layers may result in significant stress amplification in strong (elastic) layers. Finally, we emphasize the crucial role of initial stress assessment.

SummaryWe build a model of discretization errors, known as directional aliasing, to theoretically evaluate how biases in the microtremor spatial autocorrelation (SPAC) coefficient, or the real part of the SPAC spectrum of microtremor analysis, are related to the magnitudes of the imaginary part when a seismic array of only two sensors is used. By using this model, we investigate the potential utility of the imaginary spectrum component as an indicator of applicability of the two-sensor SPAC method to the field of microtremors generated at an observation site. Field data of microtremors from compact seismic arrays (1–15 m) are used to test the model. It is found that, when the imaginary components are very large in magnitude (where the threshold depends on the rk, the array radius times the wavenumber), the field of microtremors is dominated by waves arriving from a single direction parallel to the array axis and the SPAC coefficients tend to be underestimated in small rk ranges (i.e. rk < 3.8; the range considered throughout this study). In the present study, which is based on the observations of 400 microtremor arrays, the underestimates seldom exceeded 30 per cent. The SPAC coefficient estimates could be corrected in that case by using information on the imaginary part. When the imaginary components are very modest in magnitude, by contrast, there are two possible scenarios: either (i) the waves are arriving predominantly from a single direction perpendicular to the array axis and the SPAC coefficients are wildly overestimated (i.e. there was a small percentage of low-quality data, with relative errors exceeding +50 per cent, based on the observed data analyses), or (ii) the wavefield is close to isotropic and the SPAC coefficients are unbiased (i.e. 70–90 per cent of all observed data fell within the relative error range of ±20 per cent). It is difficult in that case to have SPAC coefficient estimates corrected by using information on the imaginary part alone.

SummaryAcquisition of multiple seismic datasets at different moments in time is capable of satisfying the continuously increasing demand for high-quality subsurface images to reveal both static and dynamic elements during the field development. However, in practice, challenges of pursuing this strategy lie in different perspectives related to budgetary, operational and regulatory constraints. Seismic surveys performed in a compressed manner in time and/or space can provide high-quality seismic datasets in a cost-effective and productive manner. This way of acquisition normally accompanies decompression of recorded data such as deblending and/or data reconstruction. The performance of the recovery process is of fundamental importance in determining the technical success of compressed measurements. Our proposed approach aims at realizing the benefits from compression in data acquisition, contributing to cost and efficiency, while recovering deblended and reconstructed data of sufficient quality. The approach deals jointly with deblending and data reconstruction via a sparse inversion in the frequency-wavenumber domain, coupled with constraints on causality and coherency. Additionally, we formulate a single objective function aimed at sharing static information among vintages and, at the same time, at extracting dynamic changes in the reservoir of interest according to prior subsurface information. We apply the proposed approach to both synthetic and real data. A comparison with a strategy that independently recovers compressed datasets demonstrates the viability of the proposed simultaneous method to effectively enhance the quality of recovered data and extract reliable time-lapse signatures.

SummaryKrafla is an active volcanic field and a high-temperature geothermal system in northeast Iceland. As part of a program to produce more energy from higher temperature wells, the IDDP-1 well was drilled in 2009 to reach supercritical fluid conditions below the Krafla geothermal field. However, drilling ended prematurely when the well unexpectedly encountered rhyolite magma at a depth of 2.1 km. In this paper we re-examine the magnetotelluric (MT) data that were used to model the electrical resistivity structure at Krafla. We present a new three-dimensional resistivity model that differs from previous inversions due to (1) using the full impedance tensor data and (2) a finely discretized mesh with horizontal cell dimensions of 100 m by 100 m. We obtained similar resistivity models from using two different prior models: a uniform half-space, and a previously published one-dimensional resistivity model. Our model contains a near-surface resistive layer of unaltered basalt and a low resistivity layer of hydrothermal alteration (C1). A resistive region (R1) at 1 to 2 km depth corresponds to chlorite-epidote alteration minerals that are stable at temperatures of about 220 to 500° C. A low resistivity feature (C2) coincides with the Hveragil fault system, a zone of increased permeability allowing interaction of aquifer fluids with magmatic fluids and gases. Our model contains a large, low resistivity zone (C3) below the northern half of the Krafla volcanic field that domes upward to a depth of about 1.6 km b.s.l. C3 is partially coincident with reported low S-wave velocity zones which could be due to partial melt or aqueous fluids. The low resistivity could also be attributed to dehydration and decomposition of chlorite and epidote that occurs above 500° C. As opposed to previously published resistivity models, our resistivity model shows that IDDP-1 encountered rhyolite magma near the upper edge of C3, where it intersects C2.In order to assess the sensitivity of the MT data to melt at the bottom of IDDP-1, we added hypothetical magma bodies with resistivities of 0.1 to 30 Ωm to our resistivity model and compared the synthetic MT data to the original inversion response. We used two methods to compare the MT data fit: (1) the change in r.m.s. misfit, and (2) an asymptotic p-value obtained from the Kolmogorov-Smirnov (K-S) statistical test on the two sets of data residuals. We determined that the MT data can only detect sills that are unrealistically large (2.25 km3) with very low resistivities (0.1 Ωm or 0.3 Ωm). Smaller magma bodies (0.125 km3 and 1 km3) were not detected; thus the MT data are not sensitive to small rhyolite magma bodies near the bottom of IDDP-1. Our tests gave similar results when evaluating the changes in r.m.s. misfit and the K-S test p-values, but the K-S test is a more objective method than appraising a relative change in r.m.s. misfit. Our resistivity model and resolution tests are consistent with the idea of rhyolite melt forming by re-melting of hydrothermally altered basalt on the edges of a deeper magma body.

SummaryTectonic activity is very difficult to study in the Santorini Volcanic Complex as it comprises a cluster of small/awkwardly shaped islands covered by pyroclastic deposits from which tell-tale markers are swiftly erased, while seismicity is generally absent. We address the problem by combining geophysical exploration methods to evaluate the long-term effects of tectonic deformation and time-lapse differential GPS to directly evaluate the magnitude and kinematics of present-day deformation. The former comprise 3-D gravity modelling to investigate the footprint of tectonics on the pre-volcanic Alpine basement and natural-field EM induction to map conductivity anomalies epiphenomenal to fluid circulation in faults. Our analysis identified the following principal tectonic elements:The Trans-Santorin Divide (TSD), a segmented NNW-SSE dextral strike-slip fault splitting the SVC sideways of the line joining Cape Exomytis, the Kammeni Islets and the Oia–Therassia Strait. It is collocated with a major vertical conductive zone and forms a series of dents and depressions in the basement. The Columbo Fault Zone (CFZ) is a pair of parallel NE-SW sub-vertical normal-sinistral faults straddling the northern SVC and terminating against the TSD; it may be associated with fluid injection into the shallow crust but appears to have limited effect on crustal conductivity (compared to TSD). The Anhydros Fault Zone (AFZ) is detected by its footprint on the basement, as a set of parallel northerly dipping NE-SW faults between the Athinios–Monolithos line and Fira. If it has any heave, it is left-lateral. It does not have distinguishable electrical signature and does not contribute to present-day horizontal deformation. The CFZ and AFZ are antithetic and form a graben containing the volcanic centre of Kammeni Islets.E-W extension was identified lengthwise of a zone stretching from Cape Exomytis to Athinios and along the east flank of the caldera to Imerovigli. N-S normal faulting confirmed therein, may have contributed to the localization of the east caldera wall. NNE-SSW compression was observed at SW Thera; this may have produced E-W failure and contributed to the localization of the south caldera wall. The footprint of the caldera on the basement is a parallelogram with N-S long and WNW-ESE short dimensions: if the east and south flanks collapsed along N-S normal and E-W inverse failures, then the west and north flanks may have formed analogously. Present-day deformation is localized on the TSD and CFZ: this can only be explained if the former is the synthetic (dextral) Riedel-R shear and the latter the antithetic (sinistral) Riedel-R′ shear, generated by N-S σ1 and E-W σ3 principal stress axes. Accordingly, NW-SE right-lateral shearing of the broader area is expected and indicated by several lines of indirect evidence. The geographic extent of this shearing and its role in the regional tectonics of the south Aegean remains to be confirmed and appraised by future research.Contemporary volcanic centres develop at the interface of the TSD with the CFZ/AFZ graben; volcanism appears to be controlled by tectonics and the SVC to be shaped by tectonic rather than volcanic activity.

SummaryPotential field interpretation is a powerful method to locate deep buried tectonic fault lines that contribute to intraplate earthquakes. A magnitude 5.4 earthquake (2017M5.4_PO) occurred in the Middle-Miocene Pohang Basin (PB), SE of the Korean Peninsula on 15 November 2017, in the area where no fault lines appear on geological and tectonic maps. To constrain fault locations, we calculate the gravity effect of the current basin fill with a gravity stripping method and used curvature analysis to identify former geologic and tectonic structures, assumed formed in the Early-Miocene. The Early-Miocene PB is divided into two sub-regions (northern- and southern sub-basins) by a modelled NW-SE fault line, similar to the other Early-Miocene basins (e.g. Eoil basin). Fault line trends are NE-SW in the northern sub-basin, and NNE to SSW in the southern sub-basin. 2017M5.4_PO arose from a tectonic movement along the eastern boundary of the northern sub-basin, the cross-over area from the isolated High-magnetic/Low-gravity region to Low-magnetic/High-gravity region. The largest aftershock of the 2017M5.4_PO occurred along the NW-SE fault line bordering the northern- and southern sub-basin.

SummaryNorth Korea conducted sixth underground nuclear test on September, 3rd, 2017. Unlike its previous tests, a rare subsequent collapse event occurred after about 8.5 minutes. As two types of distinctive shallow seismic events, accurate inversion of their focal mechanisms is important for event identification for CTBT. In this paper, we carry out moment tensor inversion of the nuclear test and the collapse event with gCAP using waveform data from dense regional seismic stations. And their focal mechanisms are further constrained with surface wave amplitude ratio. The results show that the surface wave amplitude ratio has further constraints for screening the waveform inversion results. The resolution of the focal mechanism inversion for the nuclear test is high, which is close to a Crack source. However, the resolution for the collapse event inversion is not so high and the source type is difficult to be accurately determined. One reason of the poor resolution for the collapse event may be due to the limited availability of high quality data, and complexity of the source process might be another factor.

SummaryCracks universally exist in Earth's crustal rocks. Many rocks are intrinsically anisotropic, which, when coupled with crack-induced anisotropy, significantly affect seismic wave propagation through the rocks. Using the method of sphere-equivalency of effective scattering, we have developed a technique to model the effective moduli of TI media containing cracks. The modeling results show the wave characteristics are significantly affected by the interaction of the two anisotropy mechanisms. To validate the validity and accuracy, the theory was applied to a recent experiment made with a VTI medium containing cracks and shows significantly better agreement with the data. For a more realistic situation, the new modeling was applied to interpret the borehole acoustic anisotropy measurement results from a fractured VTI formation, showing that the theory can adequately explain the anisotropic characteristics of the field data. With the validation and testing, the theoretical result advocated in this study can be used with confidence.

SummaryThe aim of this paper is to study the temporal variations in the seismic wavefield associated with the stress changes in the dynamic features of the Mt. Etna volcanic activity. We used Shear Wave Splitting analysis on a huge dataset of local earthquakes, in order to identify changes of the local stress field at Mt. Etna during the time interval from 2006 to 2011. This analysis allows us to obtain two parameters: the polarization direction of the fast shear wave (φ) and the time delay of the slow shear wave (Td) (time delay between the split shear waves). Orientation of φ generally provides information about the anisotropic symmetry and stress direction whereas Td provide information about the average crack density along the ray path.Based on our findings it is possible to divide Etna Volcano in three different sectors, each one distinguished by typical fast wave polarization direction. We find that the western part of the volcano is controlled by the regional tectonic stress field having a NS and EW directions. Instead, the eastern part of the volcano is mainly controlled by the local volcanic stress, particularly an EW local stress field in the NE sector (Pernicana), and a quasi NS local stress field in the SE sector (Mascalucia, Timpe), where previous studies evidenced: i) some low-Qp anomaly regions between 0 and about 6 km depth, probably associated with high pore pressure and the intense faulting (De Gori et al. 2004) and ii) by magnetotelluric surveys, several high conductivity zones, up to 8 km depth, related to a diffuse presence of hydrothermal activity and fluid circulation (Siniscalchi et al. 2012). Temporal variations in time delay, mostly before the 2008–2009 lateral eruption, can be interpreted as stress accumulation increase with a consequent release of stress due to coalescing of microcracks in the conduit for the eruption of magma.

SummarySp receiver functions have been widely used to detect the lithosphere-asthenosphere boundary (LAB) and other mantle discontinuities. However, traditional common-conversion point (CCP) stacking can be biased by the assumption of horizontal layers and this method typically underestimates scattering amplitudes from velocity boundaries with significant dips. A new pre-stack migration method based on recently developed Sp scattering kernels offers an alternative that more accurately captures the timing and amplitude of scattering. When calculating kernels, Sp-S times are estimated with the fast-marching method, and scattering amplitude versus direction, geometrical spreading, and phase shifts are accounted for. To minimize spatial aliasing with larger station spacing, Sp receiver functions are interpolated to more closely spaced pseudo-stations using either compressive sampling or spatial averaging algorithms. To test the kernel-based stacking method, synthetic Sp phases were predicted using SPECFEM2D for velocity models with a flat Moho and a negative mantle velocity gradient with a ramp structure. The kernel-based stacking method resolves horizontal interfaces equally well as CCP stacking and outperforms CCP stacking when imaging boundaries with dips of more than 8°, although dip resolution is still limited. Use of more vertically incident phases such as SKSp improves retrieval of dipping discontinuity segments. A second approach is to down-weight the portions of the kernels that have the greatest positive interference among neighboring stations, thus enhancing scattering from dipping structures where positive interference is lower. With this down-weighting, the kernel-based stacking method applied to Sp data is able to continuously resolve LAB discontinuities with dips up to 15° and to partially resolve continuous LAB discontinuities with dips of ∼20°. The intrinsic properties of teleseismic Sp phase kernels limit their ability to resolve LAB structures with dips of ∼20°-35°, but still larger dips of ∼40°-50° are resolvable with dense and appropriately placed stations. Analysis of Sp scattering kernels also explains the effectiveness of CCP stacking for quasi-horizontal interfaces.

SummaryMarchenko methods are a suite of geophysical techniques that convert seismograms of energy created by surface sources and measured by surface receivers into seismograms that would have been recorded by a virtual receiver at an arbitrary point inside the subsurface – an operation called redatuming. In principle these redatumed seismograms contain all contributions from direct, primary (singly-reflected) and multiply-reflected waves that would have been recorded by a real subsurface receiver, without requiring prior information about interfaces that generated the reflections. The potential of these methods for seismic imaging and redatuming has been demonstrated extensively in previous literature, but only using one- and two-dimensional Marchenko methods. There remain aspects of the methods that are poorly understood when applied in a three-dimensional world, so we investigate the application of Marchenko methods to three-dimensional data, subsurface structures and wavefields. We first show that for waves propagating in three dimensions, Marchenko methods can be applied to seismic data collected using both linear (so-called 2D-seismic) and areal (3D-seismic) acquisition arrays. However, for 2D acquisition arrays the Marchenko workflow requires additional dimensionality correction factors to obtain accurate solutions, even in a subsurface that only varies with depth. Without these correction factors phase errors occur in redatumed Marchenko estimates; these errors propagate through the Marchenko algorithm and create depth errors in the Marchenko images. Furthermore, applying Marchenko methods to fully three-dimensional seismic wavefields recorded by linear (2D-seismic) arrays that contain out-of-plane reflections deteriorates surface-to-subsurface Green’s function estimates with spurious energy and resulting images are less accurate than those created using ‘conventional’ imaging methods. The application of fully three-dimensional Marchenko methods using data recorded on areal arrays solves both of the above problems, creating accurately redatumed wavefields and images with reduced artifact contamination. However, it appears that source/receiver spacing at most of $\lambda _A\Big /4$ is required for accurate results using existing Marchenko methods, where λA is the dominant wavelength and in many real 3D seismic acquisition scenarios this is impractical.