ALBERT

Description: We propose a procedure for uncertainty quantification in Probabilistic Tsunami Hazard Analysis (PTHA), with a special emphasis on the uncertainty related to statistical modelling of the earthquake source in Seismic PTHA (SPTHA), and on the separate treatment of subduction and crustal earthquakes (treated as background seismicity). An event tree approach and ensemble modelling are used in spite of more classical approaches, such as the hazard integral and the logic tree. This procedure consists of four steps: (1) exploration of aleatory uncertainty through an event tree, with alternative implementations for exploring epistemic uncertainty; (2) numerical computation of tsunami generation and propagation up to a given offshore isobath; (3) (optional) site-specific quantification of inundation; (4) simultaneous quantification of aleatory and epistemic uncertainty through ensemble modelling. The proposed procedure is general and independent of the kind of tsunami source considered; however, we implement step 1, the event tree, specifically for SPTHA, focusing on seismic source uncertainty. To exemplify the procedure, we develop a case study considering seismic sources in the Ionian Sea (central-eastern Mediterranean Sea), using the coasts of Southern Italy as a target zone. The results show that an efficient and complete quantification of all the uncertainties is feasible even when treating a large number of potential sources and a large set of alternative model formulations. We also find that (i) treating separately subduction and background (crustal) earthquakes allows for optimal use of available information and for avoiding significant biases; (ii) both subduction interface and crustal faults contribute to the SPTHA, with different proportions that depend on source-target position and tsunami intensity; (iii) the proposed framework allows sensitivity and deaggregation analyses, demonstrating the applicability of the method for operational assessments.

Description: The expanding fleet of broad-band ocean-bottom seismograph (OBS) stations is facilitating the study of the structure and seismicity of oceanic plates at regional scales. For continental studies, an important tool to characterize continental crust and mantle structure is the analysis of teleseismic P receiver functions. In the oceans, however, receiver functions potentially suffer from several limiting factors that are unique to ocean sites and plate structures. In this study, we model receiver functions for a variety of oceanic lithospheric structures to investigate the possibilities and limitations of receiver functions using OBS data. Several potentially contaminating effects are examined, including pressure reverberations from the water column for various ocean-floor depths and the effects of a layer of low-velocity marine sediments. These modelling results indicate that receiver functions from OBS data are difficult to interpret in the presence of marine sediments, but shallow-water sites in subduction zone forearcs may be suitable for constraining various crustal elements around the locked megathrust fault. We propose using a complementary approach based on transfer function modelling combined with a grid search approach that bypasses receiver functions altogether and estimates model properties directly from minimally processed waveforms. Using real data examples from the Cascadia Initiative, we show how receiver and transfer functions can be used to infer seismic properties of the oceanic plate in both shallow (Cascadia forearc) and deep (Juan de Fuca Ridge) ocean settings.

Description: A new implementation of indirect boundary element method allows simulating the elastic wave propagation in complex configurations made of embedded regions that are homogeneous with irregular boundaries or flat layers. In an older implementation, each layer of a flat layered region would have been treated as a separated homogeneous region without taking into account the flat boundary information. For both types of regions, the scattered field results from fictitious sources positioned along their boundaries. For the homogeneous regions, the fictitious sources emit as in a full-space and the wave field is given by analytical Green's functions. For flat layered regions, fictitious sources emit as in an unbounded flat layered region and the wave field is given by Green's functions obtained from the discrete wavenumber (DWN) method. The new implementation allows then reducing the length of the discretized boundaries but DWN Green's functions require much more computation time than the full-space Green's functions. Several optimization steps are then implemented and commented. Validations are presented for 2-D and 3-D problems. Higher efficiency is achieved in 3-D.

Description: Though well known for layer boundaries, the use of amplitude-versus-offset (AVO) variations for non-welded boundaries like fractures is not yet investigated. Depending on the seismic wavelength used, fractures can be regarded as thin, compliant zones in rocks, in different scales. We explore the potential of multiangle AVO inversion of P-P and P-S reflections from a fracture to estimate fracture properties. We conduct laboratory experiments to measure reflection responses of dry and wet fractures. The observed P-P reflections of the wet fracture and the fracture aperture are very well predicted by the non-welded interface model. We invert the angle-dependent P-P reflectivity of the fracture to estimate both normal and tangential fracture compliances. The estimated value of the normal compliance is accurate, and it is also possible to obtain the value of the non-zero tangential compliance. We find that supplementing the information of converted P-S reflections in the AVO inversion greatly improves the estimate of the tangential compliance. The calculated compliance ratio clearly shows the existence of fluid in the fracture. This finding can be crucial for new applications in a wide range of scale—from earthquake seismology, deep and shallow seismic exploration, to non-destructive material testing.

Description: The present study evaluates the capacity of the Boom Clay as a host rock for disposal purposes, more precisely its seismic characterization, which may assess its long-term performance to store radioactive wastes. Although the formation is relatively uniform and homogeneous, there are embedded thin layers of septaria (carbonates) that may affect the integrity of the Boom Clay. Therefore, it is essential to locate these geobodies. The seismic data to characterize the Boom Clay has been acquired at the Kruibeke test site. The inversion, which allowed us to obtain the anisotropy parameters and seismic velocities of the clay, is complemented with further information such as log and laboratory data. The attenuation properties have been estimated from equivalent formations (having similar composition and seismic velocities). The inversion yields quite consistent results although the symmetry of the medium is unusual but physically possible, since the anisotropy parameter is negative. According to a time-domain calculation of the energy velocity at four frequency bands up to 900 Hz, velocity increases with frequency, a behaviour described by the Zener model. Then, we use this model to describe anisotropy and anelasticity that are implemented into the equation of motion to compute synthetic seismograms in the space–time domain. The technique is based on memory variables and the Fourier pseudospectral method. We have computed reflection coefficients of the septaria thin layer. At normal incidence, the P -wave coefficient vanishes at specific thicknesses of the layer and there is no conversion to the S wave. For example, calculations at 600 Hz show that for thicknesses of 1 m the septarium can be detected more easily since the amplitudes are higher (nearly 0.8). Converted PS waves have a high amplitude at large offsets (between 30° and 80°) and can be useful to identify the target on this basis. Moreover, we have investigated the effect of septaria embedded in the Boom Clay with several simulations, by considering a lateral partial continuity of the calcareous thin inclusions. The simulations with layers of calcareous material show continuity of the reflections even when the percentage of carbonate within the layer is very small (5–15 per cent), while for low content of the calcareous material, isolated septaria boulders generate diffraction events. We have also simulated the stacked seismic section obtained from processing of the field data. The matching between the field and synthetic sections is acceptable.

Description: A new Matched Filtering Algorithm (MFA) is proposed for detecting and analysing microseismic events recorded by downhole monitoring of hydraulic fracturing. This method requires a set of well-located template (‘parent’) events, which are obtained using conventional microseismic processing and selected on the basis of high signal-to-noise (S/N) ratio and representative spatial distribution of the recorded microseismicity. Detection and extraction of ‘child’ events are based on stacked, multichannel cross-correlation of the continuous waveform data, using the parent events as reference signals. The location of a child event relative to its parent is determined using an automated process, by rotation of the multicomponent waveforms into the ray-centred co-ordinates of the parent and maximizing the energy of the stacked amplitude envelope within a search volume around the parent's hypocentre. After correction for geometrical spreading and attenuation, the relative magnitude of the child event is obtained automatically using the ratio of stacked envelope peak with respect to its parent. Since only a small number of parent events require interactive analysis such as picking P - and S -wave arrivals, the MFA approach offers the potential for significant reduction in effort for downhole microseismic processing. Our algorithm also facilitates the analysis of single-phase child events, that is, microseismic events for which only one of the S - or P -wave arrivals is evident due to unfavourable S/N conditions. A real-data example using microseismic monitoring data from four stages of an open-hole slickwater hydraulic fracture treatment in western Canada demonstrates that a sparse set of parents (in this case, 4.6 per cent of the originally located events) yields a significant (more than fourfold increase) in the number of located events compared with the original catalogue. Moreover, analysis of the new MFA catalogue suggests that this approach leads to more robust interpretation of the induced microseismicity and novel insights into dynamic rupture processes based on the average temporal (foreshock–aftershock) relationship of child events to parents.

Description: Azimuthal anisotropy is a powerful tool to reveal information about both the present structure and past evolution of the mantle. Anisotropic images of the upper mantle are usually obtained by analysing various types of seismic observables, such as surface wave dispersion curves or waveforms, SKS splitting data, or receiver functions. These different data types sample different volumes of the earth, they are sensitive to different length scales, and hence are associated with different levels of uncertainties. They are traditionally interpreted separately, and often result in incompatible models. We present a Bayesian inversion approach to jointly invert these different data types. Seismograms for SKS and P phases are directly inverted using a cross-convolution approach, thus avoiding intermediate processing steps, such as numerical deconvolution or computation of splitting parameters. Probabilistic 1-D profiles are obtained with a transdimensional Markov chain Monte Carlo scheme, in which the number of layers, as well as the presence or absence of anisotropy in each layer, are treated as unknown parameters. In this way, seismic anisotropy is only introduced if required by the data. The algorithm is used to resolve both isotropic and anisotropic layering down to a depth of 350 km beneath two seismic stations in North America in two different tectonic settings: the stable Canadian shield (station FFC) and the tectonically active southern Basin and Range Province (station TA-214A). In both cases, the lithosphere–asthenosphere boundary is clearly visible, and marked by a change in direction of the fast axis of anisotropy. Our study confirms that azimuthal anisotropy is a powerful tool for detecting layering in the upper mantle.

Description: We have obtained new results in the statistical analysis of global earthquake catalogues with special attention to the largest earthquakes, and we examined the statistical behaviour of earthquake rate variations. These results can serve as an input for updating our recent earthquake forecast, known as the ‘Global Earthquake Activity Rate 1’ model (GEAR1), which is based on past earthquakes and geodetic strain rates. The GEAR1 forecast is expressed as the rate density of all earthquakes above magnitude 5.8 within 70 km of sea level everywhere on earth at 0.1 x 0.1 degree resolution, and it is currently being tested by the Collaboratory for Study of Earthquake Predictability. The seismic component of the present model is based on a smoothed version of the Global Centroid Moment Tensor (GCMT) catalogue from 1977 through 2013. The tectonic component is based on the Global Strain Rate Map, a ‘General Earthquake Model’ (GEM) product. The forecast was optimized to fit the GCMT data from 2005 through 2012, but it also fit well the earthquake locations from 1918 to 1976 reported in the International Seismological Centre-Global Earthquake Model (ISC-GEM) global catalogue of instrumental and pre-instrumental magnitude determinations. We have improved the recent forecast by optimizing the treatment of larger magnitudes and including a longer duration (1918–2011) ISC-GEM catalogue of large earthquakes to estimate smoothed seismicity. We revised our estimates of upper magnitude limits, described as corner magnitudes, based on the massive earthquakes since 2004 and the seismic moment conservation principle. The new corner magnitude estimates are somewhat larger than but consistent with our previous estimates. For major subduction zones we find the best estimates of corner magnitude to be in the range 8.9 to 9.6 and consistent with a uniform average of 9.35. Statistical estimates tend to grow with time as larger earthquakes occur. However, by using the moment conservation principle that equates the seismic moment rate with the tectonic moment rate inferred from geodesy and geology, we obtain a consistent estimate of the corner moment largely independent of seismic history. These evaluations confirm the above-mentioned corner magnitude value. The new estimates of corner magnitudes are important both for the forecast part based on seismicity as well as the part based on geodetic strain rates. We examine rate variations as expressed by annual earthquake numbers. Earthquakes larger than magnitude 6.5 obey the Poisson distribution. For smaller events the negative-binomial distribution fits much better because it allows for earthquake clustering.

Description: In order to improve our understanding of hazardous underground cavities, the development and collapse of a ~200 m wide salt solution mining cavity was seismically monitored in the Lorraine basin in northeastern France. The microseismic events show a swarm-like behaviour, with clustering sequences lasting from seconds to days, and distinct spatiotemporal migration. Observed microseismic signals are interpreted as the result of detachment and block breakage processes occurring at the cavity roof. Body wave amplitude patterns indicated the presence of relatively stable source mechanisms, either associated with dip-slip and/or tensile faulting. Signal overlaps during swarm activity due to short interevent times, the high-frequency geophone recordings and the limited network station coverage often limit the application of classical source analysis techniques. To overcome these shortcomings, we investigated the source mechanisms through different procedures including modelling of observed and synthetic waveforms and amplitude spectra of some well-located events, as well as modelling of peak-to-peak amplitude ratios for the majority of the detected events. We extended the latter approach to infer the average source mechanism of many swarming events at once, using multiple events recorded at a single three component station. This methodology is applied here for the first time and represents a useful tool for source studies of seismic swarms and seismicity clusters. The results obtained with different methods are consistent and indicate that the source mechanisms for at least 50 per cent of the microseismic events are remarkably stable, with a predominant thrust faulting regime with faults similarly oriented, striking NW–SE and dipping around 35°–55°. This dominance of consistent source mechanisms might be related to the presence of a preferential direction of pre-existing crack or fault structures. As an interesting byproduct, we demonstrate, for the first time directly on seismic data, that the source radiation pattern significantly controls the detection capability of a seismic station and network.

Description: We show that higher modes are an important component of high-frequency Rayleigh waves in the cross-correlations over sedimentary basins. The particle motions provide a good test for distinguishing and separating the fundamental from the first higher mode, with the fundamental mode having retrograde and the first higher mode having prograde motion in the 1–10 s period of interest. The basement depth controls the cut-off period of the first higher mode, which coincides with a rapid increase (over period) in the particle-motion ellipticity or H / V ratio of the fundamental mode. The strong higher mode we observed is not only due to the low-velocity sedimentary layer but also due to the noise sources with significant radial component such as the basin edge scattering. It is important to correctly identify the mode order when inverting the dispersion curves because misidentifying the higher mode as fundamental will lead to an anomalous high V SV velocity.

Description: We introduce a ‘double-difference’ method for the inversion for seismic wave speed structure based on adjoint tomography. Differences between seismic observations and model predictions at individual stations may arise from factors other than structural heterogeneity, such as errors in the assumed source-time function, inaccurate timings and systematic uncertainties. To alleviate the corresponding non-uniqueness in the inverse problem, we construct differential measurements between stations, thereby reducing the influence of the source signature and systematic errors. We minimize the discrepancy between observations and simulations in terms of the differential measurements made on station pairs. We show how to implement the double-difference concept in adjoint tomography, both theoretically and practically. We compare the sensitivities of absolute and differential measurements. The former provide absolute information on structure along the ray paths between stations and sources, whereas the latter explain relative (and thus higher resolution) structural variations in areas close to the stations. Whereas in conventional tomography a measurement made on a single earthquake-station pair provides very limited structural information, in double-difference tomography one earthquake can actually resolve significant details of the structure. The double-difference methodology can be incorporated into the usual adjoint tomography workflow by simply pairing up all conventional measurements; the computational cost of the necessary adjoint simulations is largely unaffected. Rather than adding to the computational burden, the inversion of double-difference measurements merely modifies the construction of the adjoint sources for data assimilation.

Description: The creation of a magnitude-homogenized catalogue is often one of the most fundamental steps in seismic hazard analysis. The process of homogenizing multiple catalogues of earthquakes into a single unified catalogue typically requires careful appraisal of available bulletins, identification of common events within multiple bulletins and the development and application of empirical models to convert from each catalogue's native scale into the required target. The database of the International Seismological Center (ISC) provides the most exhaustive compilation of records from local bulletins, in addition to its reviewed global bulletin. New open-source tools are developed that can utilize this, or any other compiled database, to explore the relations between earthquake solutions provided by different recording networks, and to build and apply empirical models in order to harmonize magnitude scales for the purpose of creating magnitude-homogeneous earthquake catalogues. These tools are described and their application illustrated in two different contexts. The first is a simple application in the Sub-Saharan Africa region where the spatial coverage and magnitude scales for different local recording networks are compared, and their relation to global magnitude scales explored. In the second application the tools are used on a global scale for the purpose of creating an extended magnitude-homogeneous global earthquake catalogue. Several existing high-quality earthquake databases, such as the ISC-GEM and the ISC Reviewed Bulletins, are harmonized into moment magnitude to form a catalogue of more than 562 840 events. This extended catalogue, while not an appropriate substitute for a locally calibrated analysis, can help in studying global patterns in seismicity and hazard, and is therefore released with the accompanying software.

Description: This paper introduces a novel approach to constructing an effective pre-conditioner for finite-difference (FD) electromagnetic modelling in geophysical applications. This approach is based on introducing an FD contraction operator, similar to one developed for integral equation formulation of Maxwell's equation. The properties of the FD contraction operator were established using an FD analogue of the energy equality for the anomalous electromagnetic field. A new pre-conditioner uses a discrete Green's function of a 1-D layered background conductivity. We also developed the formulae for an estimation of the condition number of the system of FD equations pre-conditioned with the introduced FD contraction operator. Based on this estimation, we have established that the condition number is bounded by the maximum conductivity contrast between the background conductivity and actual conductivity. When there are both resistive and conductive anomalies relative to the background, the new pre-conditioner is advantageous over using the 1-D discrete Green's function directly. In our numerical experiments with both resistive and conductive anomalies, for a land geoelectrical model with 1:10 contrast, the method accelerates convergence of an iterative method (BiCGStab) by factors of 2–2.5, and in a marine example with 1:50 contrast, by a factor of 4.6, compared to direct use of the discrete 1-D Green's function as a pre-conditioner.

Description: Reaching the global minimum of a waveform misfit function requires careful choices about the nonlinear optimization, preconditioning and regularization methods underlying an inversion. Because waveform inversion problems are susceptible to erratic convergence associated with strong nonlinearity, one or two test cases are not enough to reliably inform such decisions. We identify best practices, instead, using four seismic near-surface problems, one regional problem and two global problems. To make meaningful quantitative comparisons between methods, we carry out hundreds of inversions, varying one aspect of the implementation at a time. Comparing nonlinear optimization algorithms, we find that limited-memory BFGS provides computational savings over nonlinear conjugate gradient methods in a wide range of test cases. Comparing preconditioners, we show that a new diagonal scaling derived from the adjoint of the forward operator provides better performance than two conventional preconditioning schemes. Comparing regularization strategies, we find that projection, convolution, Tikhonov regularization and total variation regularization are effective in different contexts. Besides questions of one strategy or another, reliability and efficiency in waveform inversion depend on close numerical attention and care. Implementation details involving the line search and restart conditions have a strong effect on computational cost, regardless of the chosen nonlinear optimization algorithm.

Description: The conventional finite-difference time-domain (FDTD) method for elastic waves suffers from the staircasing error when applied to model a curved free surface because of its structured grid. In this work, an improved, stable and accurate 3-D FDTD method for elastic wave modelling on a curved free surface is developed based on the finite volume method and enlarged cell technique (ECT). To achieve a sufficiently accurate implementation, a finite volume scheme is applied to the curved free surface to remove the staircasing error; in the mean time, to achieve the same stability as the FDTD method without reducing the time step increment, the ECT is introduced to preserve the solution stability by enlarging small irregular cells into adjacent cells under the condition of conservation of force. This method is verified by several 3-D numerical examples. Results show that the method is stable at the Courant stability limit for a regular FDTD grid, and has much higher accuracy than the conventional FDTD method.

Description: Pressure solution creep (PSC) is an important elementary process in rock friction at high temperatures where solubilities of rock-forming minerals are significantly large. It significantly changes the frictional resistance and enhances time-dependent strengthening. A recent microphysical model for PSC-involved friction of clay–quartz mixtures, which can explain a transition between dilatant and non-dilatant deformation (d-nd transition), was modified here and implemented in dynamic earthquake sequence simulations. The original model resulted in essentially a kind of rate- and state-dependent friction (RSF) law, but assumed a constant friction coefficient for clay resulting in zero instantaneous rate dependency in the dilatant regime. In this study, an instantaneous rate dependency for the clay friction coefficient was introduced, consistent with experiments, resulting in a friction law suitable for earthquake sequence simulations. In addition, a term for time-dependent strengthening due to PSC was added which makes the friction law logarithmically rate-weakening in the dilatant regime. The width of the zone in which clasts overlap or, equivalently, the interface porosity involved in PSC plays a role as the state variable. Such a concrete physical meaning of the state variable is a great advantage in future modelling studies incorporating other physical processes such as hydraulic effects. Earthquake sequence simulations with different pore pressure distributions demonstrated that excess pore pressure at depth causes deeper rupture propagation with smaller slip per event and a shorter recurrence interval. The simulated ruptures were arrested a few kilometres below the point of pre-seismic peak stress at the d-nd transition and did not propagate spontaneously into the region of pre-seismic non-dilatant deformation. PSC weakens the fault against slow deformation and thus such a region cannot produce a dynamic stress drop. Dynamic rupture propagation further down to brittle-plastic transition, evidenced by geological observations, would require even smaller frictional resistance at coseismic slip rate, suggesting the importance of implementation of dynamic weakening activated at coseismic slip rates for more realistic simulation of earthquake sequences. The present models produced much smaller afterslip at deeper parts of arrested ruptures than those with logarithmic RSF laws because of a more significant rate-strengthening effect due to linearly viscous PSC. Detailed investigation of afterslip would give a clue to understand the deformation mechanism which controls shear resistance of the fault in a region of arrest of earthquake ruptures.

Description: Surface wave tomography routinely uses empirically scaled density model in the inversion of dispersion curves for shear wave speeds of the crust and uppermost mantle. An improperly selected empirical scaling relationship between density and shear wave speed can lead to unrealistic density models beneath certain tectonic formations such as sedimentary basins. Taking the Sichuan basin east to the Tibetan plateau as an example, we investigate the differences between density profiles calculated from four scaling methods and their effects on Rayleigh wave phase velocities. Analytical equations for 1-D layered models and adjoint tomography for 3-D models are used to examine the trade-off between density and S -wave velocity structures at different depth ranges. We demonstrate that shallow density structure can significantly influence phase velocities at short periods, and thereby affect the shear wave speed inversion from phase velocity data. In particular, a deviation of 25 per cent in the initial density model can introduce an error up to 5 per cent in the inverted shear velocity at middle and lower crustal depths. Therefore one must pay enough attention in choosing a proper velocity–density scaling relationship in constructing initial density model in Rayleigh wave inversion for crustal shear velocity structure.

Description: Application of the ambient noise surface wave tomography method (ANT) for determination of the upper-mantle structure requires data on long-periodic noise ( T 〉 40 s). The ANT technique implies that noise sources are distributed almost uniformly over the surface. This is practically true for short-periodic noise, however, it is not so in the case of long periods. In this paper we show that the main contribution to noise at long periods is caused by signals from earthquakes. In some cases, they may strongly distort noise cross-correlation. This leads to an incorrect determination of surface wave velocity dispersion curves. To minimize such a distortion we propose two means: (1) to use records of noise for the periods when there is no clustering of earthquakes, such as aftershocks of strong events; (2) to stack cross-correlation functions for a period of at least three years in order to achieve sufficient uniformity of earthquake locations. Validity of this approach is demonstrated by ANT results for Europe. Tomographic reconstruction of Rayleigh wave group velocities for 10–100 s measured along interstation paths was carried out in a central part of Western Europe where resolving power of the data was the highest. Locally averaged dispersion curves were inverted to vertical S -wave velocity sections in this area. The results correspond closely to known features of the structure of the region, namely: strong difference of the crust and upper-mantle structure at the opposite sides from the Tornquist–Teisseyre Line down to ~ 250 km, penetration of high-velocity material of East European Platform lithosphere under Carpathians, as well as penetration of low-velocity asthenospheric layer from the Carpathian region towards the northeast.

Description: Sensitivity analysis with synthetic models is widely used in seismic tomography as a means for assessing the spatial resolution of solutions produced by, in most cases, linear or iterative nonlinear inversion schemes. The most common type of synthetic reconstruction test is the so-called checkerboard resolution test in which the synthetic model comprises an alternating pattern of higher and lower wave speed (or some other seismic property such as attenuation) in 2-D or 3-D. Although originally introduced for application to large inverse problems for which formal resolution and covariance could not be computed, these tests have achieved popularity, even when resolution and covariance can be computed, by virtue of being simple to implement and providing rapid and intuitive insight into the reliability of the recovered model. However, checkerboard tests have a number of potential drawbacks, including (1) only providing indirect evidence of quantitative measures of reliability such as resolution and uncertainty, (2) giving a potentially misleading impression of the range of scale-lengths that can be resolved, and (3) not giving a true picture of the structural distortion or smearing that can be caused by the data coverage. The widespread use of synthetic reconstruction tests in seismic tomography is likely to continue for some time yet, so it is important to implement best practice where possible. The goal of this paper is to develop the underlying theory and carry out a series of numerical experiments in order to establish best practice and identify some common pitfalls. Based on our findings, we recommend (1) the use of a discrete spike test involving a sparse distribution of spikes, rather than the use of the conventional tightly spaced checkerboard; (2) using data coverage (e.g. ray-path geometry) inherited from the model constrained by the observations (i.e. the same forward operator or matrix), rather than the data coverage obtained by solving the forward problem through the synthetic model; (3) carrying out multiple tests using structures of different scale length; (4) taking special care with regard to what can be inferred when using synthetic structures that closely mimic what has been recovered in the observation-based model; (5) investigating the range of structural wavelengths that can be recovered using realistic levels of imposed data noise; and (6) where feasible, assessing the influence of model parametrization error, which arises from making a choice as to how structure is to be represented.

Description: We develop an approach for simulating acousto-elastic wave phenomena, including scattering from fluid–solid boundaries, where the solid is allowed to be anisotropic, with the discontinuous Galerkin method. We use a coupled first-order elastic strain-velocity, acoustic velocity–pressure formulation, and append penalty terms based on interior boundary continuity conditions to the numerical (central) flux so that the consistency condition holds for the discretized discontinuous Galerkin weak formulation. We incorporate the fluid–solid boundaries through these penalty terms and obtain a stable algorithm. Our approach avoids the diagonalization into polarized wave constituents such as in the approach based on solving elementwise Riemann problems.

Description: It is well known that large earthquakes generally trigger aftershock sequences. However, the duration of those sequences is unclear due to the gradual power-law decay with time. The triggering time is assumed to be infinite in the epidemic type aftershock sequence (ETAS) model, a widely used statistical model to describe clustering phenomena in observed earthquake catalogues. This assumption leads to the constraint that the power-law exponent p of the Omori-Utsu decay has to be larger than one to avoid supercritical conditions with accelerating seismic activity on long timescales. In contrast, seismicity models based on rate- and state-dependent friction observed in laboratory experiments predict p ≤ 1 and a finite triggering time scaling inversely to the tectonic stressing rate. To investigate this conflict, we analyse an ETAS model with finite triggering times, which allow smaller values of p . We use synthetic earthquake sequences to show that the assumption of infinite triggering times can lead to a significant bias in the maximum likelihood estimates of the ETAS parameters. Furthermore, it is shown that the triggering time can be reasonably estimated using real earthquake catalogue data, although the uncertainties are large. The analysis of real earthquake catalogues indicates mainly finite triggering times in the order of 100 days to 10 years with a weak negative correlation to the background rate, in agreement with expectations of the rate- and state-friction model. The triggering time is not the same as the apparent duration, which is the time period in which aftershocks dominate the seismicity. The apparent duration is shown to be strongly dependent on the mainshock magnitude and the level of background activity. It can be much shorter than the triggering time. Finally, we perform forward simulations to estimate the effective forecasting period, which is the time period following a mainshock, in which ETAS simulations can improve rate estimates after the occurrence of a mainshock. We find that this effective forecasting period is only in the order of 100 days for moderate mainshocks and in the order of a few years for large events, even if the underlying triggering process lasts much longer.

Description: We present a new upper-mantle tomographic model derived solely from hum seismic data. Phase correlograms between station pairs are computed to extract phase-coherent signals. Correlograms are then stacked using the time–frequency phase-weighted stack method to build-up empirical Green's functions. Group velocities and uncertainties are measured in the wide period band of 30–250 s, following a resampling approach. Less data are required to extract reliable group velocities at short periods than at long periods, and 2 yr of data are necessary to measure reliable group velocities for the entire period band. Group velocities are first regionalized and then inverted versus depth using a simulated annealing method in which the number and shape of splines that describes the S -wave velocity model are variable. The new S -wave velocity tomographic model is well correlated with models derived from earthquakes in most areas, although in India, the Dharwar craton is shallower than in other published models.

Description: Accurate prediction of ground motion intensity and its variability is an important element in seismic hazard assessment. Simulation-based ground motion prediction has become popular in light of earthquake rupture and wave-propagation modelling methods and the rapid growth in computing power. Earthquake rupture modelling needs to be physics-based and also computationally efficient. It also requires the ability to quantify the variability of finite source models for future scenario events. A generalized pseudo-dynamic source model (e.g. earthquake source statistics model) for M w 6.5–7.0 has been constructed by analysing a number of dynamic rupture models based on 1-point and 2-point statistics of kinematic source parameters. Synthetic broad-band ground motions derived from the pseudo-dynamic source model are validated against empirical ground motion prediction equations (GMPEs). The pseudo-dynamic source model produces ground motion intensities mostly compatible with the empirical GMPEs in the broad frequency range. The perturbation analysis of 1-point and 2-point statistics helps to elucidate the effect of source statistics on ground motions in a systematic way. The constructed pseudo-dynamic source model may be used to generate a feasible range of rupture scenarios for future events and to simulate expected ground motions for seismic hazard assessment.

Description: We evaluate the bias in parameter estimates of the ETAS model. We show that when a simulated catalogue is magnitude-truncated there is considerable bias, whereas when it is not truncated there is no discernible bias. We also discuss two further implied assumptions in the ETAS and other self-exciting models. First, that the triggering boundary magnitude is equivalent to the catalogue completeness magnitude. Secondly, the assumption in the Gutenberg–Richter relationship that numbers of events increase exponentially as magnitude decreases. These two assumptions are confounded with the magnitude truncation effect. We discuss the effect of these problems on analyses of real earthquake catalogues.

Description: The frequency dependence of the quality factor Q has long been predicted by mathematical modelling and laboratory measurements; however, in situ evidence from seismic surveys is still lacking. We have conducted the cross-hole seismic surveys to investigate the near-surface seismic attenuation in the Daqing oilfield in northeastern China. The seismic waves were fired in a source hole of 40 m from the bottom to the surface at an interval of 1 m and were recorded in a receiver hole of 40 m by two geophones with one at the surface and the other one at the bottom. The direct waves were extracted to avoid the noise disturbance and the reflection interference, and the attenuations without the effects of the source signature and the receiver coupling were estimated by a method we proposed. The nonlinear attenuations were observed and fitted using the power-law-based Q . The reliability of Q estimate was verified by the high similarity between the real and the simulated attenuations. Therefore, the experiment we have conducted can be treated as a reliable evidence for the frequency dependence of near-surface  Q .

Description: The southeastern European cities of Sofia and Thessaloniki are explored as example site-specific scenarios by geographically zoning their individual localized seismic sources based on the highest probabilities of magnitude exceedance. This is with the aim of determining the major components contributing to each city's seismic hazard. Discrete contributions from the selected input earthquake catalogue are investigated to determine those areas that dominate each city's prevailing seismic hazard with respect to magnitude and source-to-site distance. This work is based on an earthquake catalogue developed and described in a previously published paper by the author and components of a magnitude probability density function. Binned magnitude and distance classes are defined using a joint magnitude–distance distribution. The prevailing seismicity to each city—as defined by a child data set extracted from the parent earthquake catalogue for each city considered—is divided into distinct constrained data bins of small discrete magnitude and source-to-site distance intervals. These are then used to describe seismic hazard in terms of uni-variate modal values; that is, M * and D * which are the modal magnitude and modal source-to-site distance in each city's local historical seismicity. This work highlights that Sofia's dominating seismic hazard—that is, the modal magnitudes possessing the highest probabilities of occurrence—is located in zones confined to two regions at 60–80 km and 170–180 km from this city, for magnitude intervals of 5.75–6.00 M w and 6.00–6.25 M w respectively. Similarly, Thessaloniki appears prone to highest levels of hazard over a wider epicentral distance interval, from 80 to 200 km in the moment magnitude range 6.00–6.25 M w .

Description: We collect two months of ambient noise data recorded by 35 broad-band seismic stations in a 9 x 11 km area (1–3 km station interval) near Karamay, China, and do cross-correlation of noise data between all station pairs. Array beamforming analysis of the ambient noise data shows that ambient noise sources are unevenly distributed and the most energetic ambient noise mainly comes from azimuths of 40°–70°. As a consequence of the strong directional noise sources, surface wave components of the cross-correlations at 1–5 Hz show clearly azimuthal dependence, and direct dispersion measurements from cross-correlations are strongly biased by the dominant noise energy. This bias renders that the dispersion measurements from cross-correlations do not accurately reflect the interstation velocities of surface waves propagating directly from one station to the other, that is, the cross-correlation functions do not retrieve empirical Green's functions accurately. To correct the bias caused by unevenly distributed noise sources, we adopt an iterative inversion procedure. The iterative inversion procedure, based on plane-wave modeling, includes three steps: (1) surface wave tomography, (2) estimation of ambient noise energy and biases and (3) phase velocities correction. First, we use synthesized data to test the efficiency and stability of the iterative procedure for both homogeneous and heterogeneous media. The testing results show that: (1) the amplitudes of phase velocity bias caused by directional noise sources are significant, reaching ~2 and ~10 per cent for homogeneous and heterogeneous media, respectively; (2) phase velocity bias can be corrected by the iterative inversion procedure and the convergence of inversion depends on the starting phase velocity map and the complexity of the media. By applying the iterative approach to the real data in Karamay, we further show that phase velocity maps converge after 10 iterations and the phase velocity maps obtained using corrected interstation dispersion measurements are more consistent with results from geology surveys than those based on uncorrected data. As ambient noise in high-frequency band (〉1 Hz) is mostly related to human activities or climate events, both of which have strong directivity, the iterative approach demonstrated here helps improve the accuracy and resolution of ANT in imaging shallow earth structures.

Description: We developed an improved method for the separation of intrinsic and scattering attenuation of seismic shear waves by envelope inversion called Qopen . The method optimizes the fit between Green's functions for the acoustic, isotropic radiative transfer theory and observed energy densities of earthquakes. The inversion allows the determination of scattering and intrinsic attenuation, site corrections and spectral source energies for the investigated frequency bands. Source displacement spectrum and the seismic moment of the analysed events can be estimated from the obtained spectral source energies. We report intrinsic and scattering attenuation coefficients of shear waves near three geothermal reservoirs in Germany for frequencies between 1 and 70 Hz. The geothermal reservoirs are located in Insheim, Landau (both Upper Rhine Graben) and Unterhaching (Molasse basin). We compare these three sedimentary sites to two sites located in crystalline rock with respect to scattering and intrinsic attenuation. The inverse quality factor for intrinsic attenuation is constant in sediments for frequencies smaller than $10\hspace{1.00006pt}\mathrm{Hz}$ and decreasing for higher frequencies. For crystalline rock, it is on a lower level and strictly monotonic decreasing with frequency. Intrinsic attenuation dominates scattering except for the Upper Rhine Graben, where scattering is dominant for frequencies below $10\hspace{1.00006pt}\mathrm{Hz}$ . Observed source displacement spectra show a high-frequency fall-off greater than or equal to 3.

Description: This paper presents a generalized wave equation which unifies viscoelastic and pure elastic cases into a single wave equation. In the generalized wave equation, the degree of viscoelasticity varies between zero and unity, and is defined by a controlling parameter. When this viscoelastic controlling parameter equals to 0, the viscous property vanishes and the generalized wave equation becomes a pure elastic wave equation. When this viscoelastic controlling parameter equals to 1, it is the Stokes equation made up of a stack of pure elastic and Newtonian viscous models. Given this generalized wave equation, an analytical solution is derived explicitly in terms of the attenuation and the velocity dispersion. It is proved that, for any given value of the viscoelastic controlling parameter, the attenuation component of this generalized wave equation perfectly satisfies the power laws of frequency. Since the power laws are the fundamental characteristics in physical observations, this generalized wave equation can well represent seismic wave propagation through subsurface media.

Description: We investigate the relationship between subduction processes and related seismicity for the Lesser Antilles Arc using the Gutenberg–Richter law. This power law describes the earthquake-magnitude distribution, with the gradient of the cumulative magnitude distribution being commonly known as the b -value. The Lesser Antilles Arc was chosen because of its along-strike variability in sediment subduction and the transition from subduction to strike-slip movement towards its northern and southern ends. The data are derived from the seismicity catalogues from the Seismic Research Centre of The University of the West Indies and the Observatoires Volcanologiques et Sismologiques of the Institut de Physique du Globe de Paris and consist of subcrustal events primarily from the slab interface. The b -value is found using a Kolmogorov–Smirnov test for a maximum-likelihood straight line-fitting routine. We investigate spatial variations in b -values using a grid-search with circular cells as well as an along-arc projection. Tests with different algorithms and the two independent earthquake cataloges provide confidence in the robustness of our results. We observe a strong spatial variability of the b -value that cannot be explained by the uncertainties. Rather than obtaining a simple north–south b -value distribution suggestive of the dominant control on earthquake triggering being water released from the sedimentary cover on the incoming American Plates, or a b -value distribution that correlates with on the obliquity of subduction, we obtain a series of discrete, high b -value ‘bull's-eyes’ along strike. These bull's-eyes, which indicate stress release through a higher fraction of small earthquakes, coincide with the locations of known incoming oceanic fracture zones on the American Plates. We interpret the results in terms of water being delivered to the Lesser Antilles subduction zone in the vicinity of fracture zones providing lubrication and thus changing the character of the related seismicity. Our results suggest serpentinization around mid-ocean ridge transform faults, which go on to become fracture zones on the incoming plate, plays a significant role in the delivery of water into the mantle at subduction zones.

Description: Given a moment tensor m inferred from seismic data for an earthquake, we define ${\scr P}(V)$ to be the probability that the true moment tensor for the earthquake lies in the neighbourhood of m that has fractional volume V . The average value of ${\scr P}(V)$ is then a measure of our confidence in  m . The calculation of ${\scr P}(V)$ requires knowing both the probability $\skew4\hat{P}(\omega )$ and the fractional volume $\skew4\hat{V}(\omega )$ of the set of moment tensors within a given angular radius of  m . We explain how to construct $\skew4\hat{P}(\omega )$ from a misfit function derived from seismic data, and we show how to calculate $\skew4\hat{V}(\omega )$ , which depends on the set $\mathbb {M}$ of moment tensors under consideration. The two most important instances of $\mathbb {M}$ are where $\mathbb {M}$ is the set of all moment tensors of fixed norm, and where $\mathbb {M}$ is the set of all double couples of fixed norm.

Description: The M s ~ 7.7 Sarez-Pamir earthquake of 1911 February 18 is the largest instrumentally recorded earthquake in the Pamir region. It triggered one of the largest landslides of the past century, building a giant natural dam and forming Lake Sarez. As for many strong earthquakes from that time, information about source parameters of the Sarez-Pamir earthquake is limited due to the sparse observations. Here, we present the analysis of analogue seismic records of the Sarez-Pamir earthquake. We have collected, scanned and digitized 26 seismic records from 13 stations worldwide to relocate the epicentre and determine the event's depth (~26 km) and magnitude ( m B 7.3 and M s 7.7). The unusually good quality of the digitized waveforms allowed their modelling, revealing an NE-striking sinistral strike-slip focal mechanism in accordance with regional tectonics. The shallow depth and magnitude ( M w 7.3) of the earthquake were confirmed. Additionally, we investigated the possible contribution of the landslide to the waveforms and present an alternative source model assuming the landslide and earthquake occurred in close sequence.

Description: To obtain the synthetic seismogram using the Cagniard-de Hoop method, one needs to calculate the integral over slowness. When the source is shallow and the slowness is near the zero of the Rayleigh function, the integrand behaves like a sharp pulse. In this study, we attempt to study this pulse with an asymptotic approach, and conclude that the Rayleigh wave in the time domain originates from this pulse in the slowness domain. We therefore offer an explanation of the excitation of the Rayleigh wave in a mathematical point of view. In addition, we propose a method to improve the efficiency of the numerical quadrature in the calculation of the synthetic seismogram.

Description: Many studies have sought to seismically image plumes rising from the deep mantle in order to settle the debate about their presence and role in mantle dynamics, yet the predicted seismic signature of realistic plumes remains poorly understood. By combining numerical simulations of flow, mineral-physics constraints on the relationships between thermal anomalies and wave speeds, and spectral-element method based computations of seismograms, we estimate the delay times of teleseismic S and P waves caused by thermal plumes. Wave front healing is incomplete for seismic periods ranging from 10 s (relevant in traveltime tomography) to 40 s (relevant in waveform tomography). We estimate P -wave delays to be immeasurably small (〈0.3 s). S -wave delays are larger than 0.4 s even for S waves crossing the conduits of the thinnest thermal plumes in our geodynamic models. At longer periods (〉20 s), measurements of instantaneous phase misfit may be more useful in resolving narrow plume conduits. To detect S -wave delays of 0.4–0.8 s and the diagnostic frequency dependence imparted by plumes, it is key to minimize the influence of the heterogeneous crust and upper mantle. We argue that seismic imaging of plumes will advance significantly if data from wide-aperture ocean-bottom networks were available since, compared to continents, the oceanic crust and upper mantle are relatively simple.

Description: Field experiments are used to unequivocally demonstrate seismic superresolution imaging of subwavelength objects in the near-field region of the source. The field test is for a conventional hammer source striking a metal plate near subwavelength scatterers and the seismic data are recorded by vertical-component geophones in the far-field locations of the sources. Time-reversal mirrors (TRMs) are then used to refocus the scattered energy with subwavelength resolution to the position of the original source. A spatial resolution of /10, where is the dominant wavelength associated with the data, is seen in the field tests that exceeds the Abbe resolution limit of /2.

Description: We propose an adaptive root-determining strategy that is very useful when dealing with trapped modes or Stoneley modes whose energies become very insignificant on the free surface in the presence of low-velocity layers or fluid layers in the model. Loss of modes in these cases or inaccuracy in the calculation of these modes may then be easily avoided. Built upon the generalized reflection/transmission coefficients, the concept of ‘family of secular functions’ that we herein call ‘adaptive mode observers’ is thus naturally introduced to implement this strategy, the underlying idea of which has been distinctly noted for the first time and may be generalized to other applications such as free oscillations or applied to other methods in use when these cases are encountered. Additionally, we have made further improvements upon the generalized reflection/transmission coefficient method; mode observers associated with only the free surface and low-velocity layers (and the fluid/solid interface if the model contains fluid layers) are adequate to guarantee no loss and high precision at the same time of any physically existent modes without excessive calculations. Finally, the conventional definition of the fundamental mode is reconsidered, which is entailed in the cases under study. Some computational aspects are remarked on. With the additional help afforded by our superior root-searching scheme and the possibility of speeding calculation using a less number of layers aided by the concept of ‘turning point’, our algorithm is remarkably efficient as well as stable and accurate and can be used as a powerful tool for widely related applications.

Description: Numerical solvers of wave equations have been widely used to simulate global seismic waves including PP waves for modelling 410/660 km discontinuity and Rayleigh waves for imaging crustal structure. In order to avoid extra computation cost due to ocean water effects, these numerical solvers usually adopt water column approximation, whose accuracy depends on frequency and needs to be investigated quantitatively. In this paper, we describe a unified representation of accurate and approximate forms of the equivalent water column boundary condition as well as the free boundary condition. Then we derive an analytical form of the PP -wave reflection coefficient with the unified boundary condition, and quantify the effects of water column approximation on amplitude and phase shift of the PP waves. We also study the effects of water column approximation on phase velocity dispersion of the fundamental mode Rayleigh wave with a propagation matrix method. We find that with the water column approximation: (1) The error of PP amplitude and phase shift is less than 5 per cent and 9° at periods greater than 25 s for most oceanic regions. But at periods of 15 s or less, PP is inaccurate up to 10 per cent in amplitude and a few seconds in time shift for deep oceans. (2) The error in Rayleigh wave phase velocity is less than 1 per cent at periods greater than 30 s in most oceanic regions, but the error is up to 2 per cent for deep oceans at periods of 20 s or less. This study confirms that the water column approximation is only accurate at long periods and it needs to be improved at shorter periods.

Description: Seismic wave resonance in sedimentary basins is a well-recognized seismic hazard; however, concentrated areas of earthquake damage have been observed near basin edges, where wave propagation is particularly complex and difficult to understand with sparse observations. The Tokyo metropolitan area is densely populated, subject to strong shaking from a diversity of earthquake sources, and sits atop the deep Kanto sedimentary basin. It is also instrumented with two seismic arrays: the dense MEtropolitan Seismic Observation network (MeSO-net) within the basin, and the High sensitivity seismograph network (Hi-net) surrounding it. In this study, we explore the 3-D seismic wavefield within and throughout the Kanto basin, including near and across basin boundaries, using cross-correlations of all components of ambient seismic field between the stations of these two arrays. Dense observations allow us to observe clearly the propagation of three modes of both Rayleigh and Love waves. They also show how the wavefield behaves in the vicinity of sharp basin edges with reflected/converted waves and excitation of higher modes.

Description: We present an analytical approach to compute the curvature effect by the new analytical solutions of coseismic deformation derived for the homogeneous sphere model. We consider two spheres with different radii: one is the same as earth and the other with a larger radius can approximate a half-space model. Then, we calculate the coseismic displacements for the two spheres and define the relative percentage of the displacements as the curvature effect. The near-field curvature effect is defined relative to the maximum coseismic displacement. The results show that the maximum curvature effect is about 4 per cent for source depths of less than 100 km, and about 30 per cent for source depths of less than 600 km. For the far-field curvature effect, we define it relative to the observing point. The curvature effect is extremely large and sometimes exceeds 100 per cent. Moreover, this new approach can be used to estimate any planet's curvature effect quantitatively. For a smaller sphere, such as the Moon, the curvature effect is much larger than that of the Earth, with an inverse ratio to the earth's radius.

Description: In Antarctica, locally grounded ice, such as ice rises bordering floating ice shelves, plays a major role in the ice mass balance as it stabilizes the ice sheet flow from the hinterland. When in direct contact with the ocean, the ice rise buttressing effect may be altered in response of changing ocean forcing. To investigate this vulnerable zone, four sites near the boundary of an ice shelf with an ice rise promontory in Dronning Maud Land, East-Antarctica were monitored for a month in early 2014 with new instruments that include both seismic and GPS sensors. Our study indicated that this transition zone experiences periodic seismic activity resulting from surface crevassing during oceanic tide-induced flexure of the ice shelf. The most significant finding is the observation of apparent fortnightly tide-modulated low-frequency, long-duration seismic events at the seaward front of the ice rise promontory. A basal origin of these events is postulated with the ocean water surge at each new spring tide triggering basal crevassing or basal slip on a local bedrock asperity. Detection and monitoring of such seismicity may help identifying ice rise zones vulnerable to intensified ocean forcing.

Description: SKS arrivals from ocean bottom seismometer (OBS) data from an offshore southern California deployment are analysed for shear wave splitting. The project involved 34 OBSs deployed for 12 months in a region extending up to 500 km west of the coastline into the oceanic Pacific plate. The measurement process consisted of removing the effects of anisotropy using a range of values for splitting fast directions and delay times to minimize energy along the transverse seismometer axis. Computed splitting parameters are unexpectedly similar to onland parameters, exhibiting WSW–ENE fast polarization directions and delays between 0.8 and 1.8 s, even for oceanic plate sites. This is the first SKS splitting study to extend across the entire boundary between the North America and Pacific plates, into the oceanic part of the Pacific plate. The splitting results show that the fast direction of anisotropy on the Pacific plate does not align with absolute plate motion (APM), and they extend the trend of anisotropy in southern California an additional 500 km west, well onto the oceanic Pacific plate. We model the finite strain and anisotropy within the asthenosphere associated with density–buoyancy driven mantle flow and the effects of APM. In the absence of plate motion effects, such buoyancy driven mantle flow would be NE-directed beneath the Pacific plate observations. The best-fit patterns of mantle flow are inferred from the tomography-based models that show primary influences from foundering higher-density zones associated with the history of subduction beneath North America. The new offshore SKS measurements, when combined with measurements onshore within the plate boundary zone, indicate that dramatic lateral variations in density-driven upper-mantle flow are required from offshore California into the plate boundary zone in California and western Basin and Range.

Description: The islands composing the Maltese archipelago (Central Mediterranean) are characterized by a four-layer sequence of limestones and clays. A common feature found in the western half of the archipelago is Upper Coralline Limestone (UCL) plateaus and hillcaps covering a soft Blue Clay (BC) layer which can be up to 75 m thick. The BC layer introduces a velocity inversion in the stratigraphy, implying that the V S 30 (traveltime average sear wave velocity ( V S ) in the upper 30 m) parameter is not always suitable for seismic microzonation purposes. Such a layer may produce amplification effects, however might not be included in the V S 30 calculations. In this investigation, V S profiles at seven sites characterized by such a lithological sequence are obtained by a joint inversion of the single-station Horizontal-to-Vertical Spectral Ratios (H/V or HVSR) and effective dispersion curves from array measurements analysed using the Extended Spatial Auto-Correlation technique. The lithological sequence gives rise to a ubiquitous H/V peak between 1 and 2 Hz. All the effective dispersion curves obtained exhibit a ‘normal’ dispersive trend at low frequencies, followed by an inverse dispersive trend at higher frequencies. This shape is tentatively explained in terms of the presence of higher mode Rayleigh waves, which are commonly present in such scenarios. Comparisons made with the results obtained at the only site in Malta where the BC is missing below the UCL suggest that the characteristics observed at the other seven sites are due to the presence of the soft layer. The final profiles reveal a variation in the V S of the clay layer with respect to the depth of burial and some regional variations in the UCL layer. This study presents a step towards a holistic seismic risk assessment that includes the implications on the site effects induced by the buried clay layer. Such assessments have not yet been done for Malta.

Description: Noise is a persistent feature in seismic data and so poses challenges in extracting increased accuracy in seismic images and physical interpretation of the subsurface. In this paper, we analyse passive seismic data from the Aquistore carbon capture and storage pilot project permanent seismic array to characterise, classify and model seismic noise. We perform noise analysis for a three-month subset of passive seismic data from the array and provide conclusive evidence that the noise field is not white, stationary, or Gaussian; characteristics commonly yet erroneously assumed in most conventional noise models. We introduce a novel noise modelling method that provides a significantly more accurate characterisation of real seismic noise compared to conventional methods, which is quantified using the Mann–Whitney–White statistical test. This method is based on a statistical covariance modelling approach created through the modelling of individual noise signals. The identification of individual noise signals, broadly classified as stationary, pseudo-stationary and non-stationary, provides a basis on which to build an appropriate spatial and temporal noise field model. Furthermore, we have developed a workflow to incorporate realistic noise models within synthetic seismic data sets providing an opportunity to test and analyse detection and imaging algorithms under realistic noise conditions.

Description: Plate-scale deformation is expected to impart seismic anisotropic fabrics on the lithosphere. Determination of the fast shear wave orientation ( ) and the delay time between the fast and slow split shear waves ( t ) via SKS splitting can help place spatial and temporal constraints on lithospheric deformation. The Canadian Appalachians experienced multiple episodes of deformation during the Phanerozoic: accretionary collisions during the Palaeozoic prior to the collision between Laurentia and Gondwana, and rifting related to the Mesozoic opening of the North Atlantic. However, the extent to which extensional events have overprinted older orogenic trends is uncertain. We address this issue through measurements of seismic anisotropy beneath the Canadian Appalachians, computing shear wave splitting parameters ( , t ) for new and existing seismic stations in Nova Scotia and New Brunswick. Average t values of 1.2 s, relatively short length scale (≥100 km) splitting parameter variations, and a lack of correlation with absolute plate motion direction and mantle flow models, demonstrate that fossil lithospheric anisotropic fabrics dominate our results. Most fast directions parallel Appalachian orogenic trends observed at the surface, while t values point towards coherent deformation of the crust and mantle lithosphere. Mesozoic rifting had minimal impact on our study area, except locally within the Bay of Fundy and in southern Nova Scotia, where fast directions are subparallel to the opening direction of Mesozoic rifting; associated t values of 〉1 s require an anisotropic layer that spans both the crust and mantle, meaning the formation of the Bay of Fundy was not merely a thin-skinned tectonic event.

Description: A new method for discrete sampling of signals is presented with specific applications to the reconstruction of recorded interfering wavefields from two or more sources excited simultaneously at discrete positions along lines. By utilizing a periodic sequence of source signatures along one of the source lines, the corresponding wavefield becomes separately visible in a part of the spectral domain where it can be isolated, processed and subtracted from the interfering wavefields. As a result, interfering wavefields from multiple sources recorded at a single location can be fully separated from each other. The concept is referred to as signal apparition which we suggest refers to ‘the act of becoming visible’. It may find applications in a wide range of disciplines relying on wave experimentation, such as acoustic, seismic and electromagnetic imaging of the Earth's interior for instance to significantly enhance resolution of subsurface images.

Description: In the context of seismic monitoring, recent studies made successful use of seismic coda waves to locate medium changes on the horizontal plane. Locating the depth of the changes, however, remains a challenge. In this paper, we use 3-D wavefield simulations to address two problems: first, we evaluate the contribution of surface- and body-wave sensitivity to a change at depth. We introduce a thin layer with a perturbed velocity at different depths and measure the apparent relative velocity changes due to this layer at different times in the coda and for different degrees of heterogeneity of the model. We show that the depth sensitivity can be modelled as a linear combination of body- and surface-wave sensitivity. The lapse-time-dependent sensitivity ratio of body waves and surface waves can be used to build 3-D sensitivity kernels for imaging purposes. Second, we compare the lapse-time behaviour in the presence of a perturbation in horizontal and vertical slabs to address, for instance, the origin of the velocity changes detected after large earthquakes.

Description: Traveltime tomography is the main method by which the Earth's seismic velocity is determined on all scales, from the near-surface (〈100 m) to the core. Usually traveltime tomography uses ray theory, an infinite-frequency approximation of wave propagation. A theory developed in global seismology to account for the finite-frequency nature of seismic data, known as finite-frequency traveltime tomography (FFTT), can theoretically provide a more accurate estimation of velocity. But the FFTT theory is generally not applicable to near-surface data because there is no reference velocity model known in advance that is capable of yielding synthetic waveforms that are close enough to the recorded seismograms to yield a reliable delay time. Also, there is usually no reference model for which the unknown velocity model represents a small (linear) perturbation from the reference model. This paper presents a frequency dependent form of non-linear traveltime tomography specifically designed for near-surface seismic data in which a starting model, iterative approach with recalculated travel paths at each iteration, and the calculation of a frequency-dependent total traveltime, as opposed to a delay time, are used. Frequency-dependent traveltime tomography (FDTT) involves two modifications to conventional traveltime tomography: (1) the calculation of frequency-dependent traveltimes using wavelength-dependent velocity smoothing (WDVS) and (2) the corresponding sensitivity kernels that arise from using WDVS. Results show that the former modification is essential to achieve significant benefits from FDTT, whereas the latter is optional in that similar results can be achieved using infinite-frequency kernels. The long seismic wavelengths relative to the total path lengths and the size of subsurface heterogeneities of typical near-surface data means the improvements over ray theory tomography are significant. The benefits of FDTT are demonstrated using conventional minimum-structure regularization techniques to address the issue of model non-uniqueness. For synthetic data, the estimated FDTT models are shown to be more accurate than the corresponding infinite-frequency-derived models. Both 2-D and 3-D applications of FDTT to real data from a near-surface study yield estimated models that contain more structure than the corresponding infinite-frequency-derived models. Applications of FDTT without regularization demonstrate the potential of the WDVS-derived sensitivity kernels to provide a natural smoothing of the velocity model and thereby allow the data alone to determine the final model structure.

Description: We analysed a catalogue of Italian earthquakes, covering 55 yr of data from 1960 to 2014 with magnitudes homogeneously converted to M w , to compute the time-dependent relative frequencies with which strong seismic shocks (4.0 ≤ M w 〈 5.0), widely felt by the population, have been followed by main shocks ( M w ≥ 5.0) that threatened the health and the properties of the persons living in the epicentral area. Assuming the stationarity of the seismic release properties, such frequencies are estimates of the probabilities of potentially destructive shocks after the occurrence of future strong shocks. We compared them with the time-independent relative frequencies of random occurrence in terms of the frequency gain that is the ratio between the time-dependent and time-independent relative frequencies. The time-dependent relative frequencies vary from less than 1 per cent to about 20 per cent, depending on the magnitudes of the shocks and the time windows considered (ranging from minutes to years). They remain almost constant for a few hours after the strong shock and then decrease with time logarithmically. Strong earthquakes (with M w ≥ 6.0) mainly occurred within two or three months of the strong shock. The frequency gains vary from about 10 000 for very short time intervals to less than 10 for a time interval of 2 yr. Only about 1/3 of main shocks were preceded by at least a strong shock in the previous day and about 1/2 in the previous month.

Description: We present the first high-resolution Rayleigh-wave phase-velocity azimuthal anisotropy tomography of the Japan subduction zone at periods of 20–150 s, which is determined using a large number of high-quality amplitude and phase data of teleseismic fundamental-mode Rayleigh waves. The obtained 2-D anisotropic phase-velocity models are then inverted for a 3-D shear-wave velocity azimuthal anisotropy tomography down to a depth of ~300 km beneath Japan. The subducting Pacific slab is imaged as a dipping high-velocity zone with trench-parallel fast-velocity directions (FVDs) which may indicate the anisotropy arising from the normal faults produced at the outer-rise area near the Japan trench axis, overprinting the slab fossil fabric, whereas the mantle wedge generally exhibits lower velocities with trench-normal FVDs which reflect subduction-driven corner flow and anisotropy. Depth variations of azimuthal anisotropy are revealed in the big mantle wedge beneath the Japan Sea, which may reflect past deformations in the Eurasian lithosphere related to backarc spreading during 21 to 15 Ma and complex current convection in the asthenosphere induced by active subductions of both the Pacific and Philippine Sea plates.

Description: Markov chain Monte Carlo sampling methods are widely used for non-linear Bayesian inversion where no analytical expression for the forward relation between data and model parameters is available. Contrary to the linear(ized) approaches, they naturally allow to evaluate the uncertainties on the model found. Nevertheless their use is problematic in high-dimensional model spaces especially when the computational cost of the forward problem is significant and/or the a posteriori distribution is multimodal. In this case, the chain can stay stuck in one of the modes and hence not provide an exhaustive sampling of the distribution of interest. We present here a still relatively unknown algorithm that allows interaction between several Markov chains at different temperatures. These interactions (based on importance resampling) ensure a robust sampling of any posterior distribution and thus provide a way to efficiently tackle complex fully non-linear inverse problems. The algorithm is easy to implement and is well adapted to run on parallel supercomputers. In this paper, the algorithm is first introduced and applied to a synthetic multimodal distribution in order to demonstrate its robustness and efficiency compared to a simulated annealing method. It is then applied in the framework of first arrival traveltime seismic tomography on real data recorded in the context of hydraulic fracturing. To carry out this study a wavelet-based adaptive model parametrization has been used. This allows to integrate the a priori information provided by sonic logs and to reduce optimally the dimension of the problem.

Description: Using data from more than 2000 seismic stations from multiple networks arrayed throughout China (CEArray, China Array, NECESS, PASSCAL, GSN) and surrounding regions (Korean Seismic Network, F-Net, KNET), we perform ambient noise Rayleigh wave tomography across the entire region and earthquake tomography across parts of South China and Northeast China. We produce isotropic Rayleigh wave group and phase speed maps with uncertainty estimates from 8 to 50 s period across the entire region of study, and extend them to 70 s period where earthquake tomography is performed. Maps of azimuthal anisotropy are estimated simultaneously to minimize anisotropic bias in the isotropic maps, but are not discussed here. The 3D model is produced using a Bayesian Monte Carlo formalism covering all of China, extending eastwards through the Korean Peninsula, into the marginal seas, to Japan. We define the final model as the mean and standard deviation of the posterior distribution at each location on a 0.5° x 0.5° grid from the surface to 150 km depth. Surface wave dispersion data do not strongly constrain internal interfaces, but shear wave speeds between the discontinuities in the crystalline crust and uppermost mantle are well determined. We design the resulting model as a reference model, which is intended to be useful to other researchers as a starting model, to predict seismic wave fields and observables and to predict other types of data (e.g. topography, gravity). The model and the data on which it is based are available for download. In addition, the model displays a great variety and considerable richness of geological and tectonic features in the crust and in the uppermost mantle deserving of further focus and continued interpretation.

Description: Earthquake source inversion is highly dependent on location determination and velocity models. Uncertainties in both the model parameters and the observations need to be rigorously incorporated into an inversion approach. Here, we show a probabilistic Bayesian method that allows formal inclusion of the uncertainties in the moment tensor inversion. This method allows the combination of different sets of far-field observations, such as P -wave and S -wave polarities and amplitude ratios, into one inversion. Additional observations can be included by deriving a suitable likelihood function from the uncertainties. This inversion produces samples from the source posterior probability distribution, including a best-fitting solution for the source mechanism and associated probability. The inversion can be constrained to the double-couple space or allowed to explore the gamut of moment tensor solutions, allowing volumetric and other non-double-couple components. The posterior probability of the double-couple and full moment tensor source models can be evaluated from the Bayesian evidence, using samples from the likelihood distributions for the two source models, producing an estimate of whether or not a source is double-couple. Such an approach is ideally suited to microseismic studies where there are many sources of uncertainty and it is often difficult to produce reliability estimates of the source mechanism, although this can be true of many other cases. Using full-waveform synthetic seismograms, we also show the effects of noise, location, network distribution and velocity model uncertainty on the source probability density function. The noise has the largest effect on the results, especially as it can affect other parts of the event processing. This uncertainty can lead to erroneous non-double-couple source probability distributions, even when no other uncertainties exist. Although including amplitude ratios can improve the constraint on the source probability distribution, the measurements are often systematically affected by noise, leading to deviation from their noise-free true values and consequently adversely affecting the source probability distribution, especially for the full moment tensor model. As an example of the application of this method, four events from the Krafla volcano in Iceland are inverted, which show clear differentiation between non-double-couple and double-couple sources, reflected in the posterior probability distributions for the source models.

Description: Extracting accurate empirical Green's functions from the ambient seismic noise field requires the noise to be fully diffuse and that different frequency components are not correlated. Calculating a matrix of correlation coefficients of power spectral samples can be used to estimate deviations from a fully diffuse random noise field in the analysed frequency range. A fully diffuse field has correlations only in a narrow region around the diagonal of the matrix, with frequency resolution inversely proportional to length of the used time window. Analysis of low-frequency data (0.005–0.6 Hz) recorded by three broad-band stations of the southern California seismic network reveals three common types of correlations, manifested in the correlation coefficient matrix as square, diagonal halo and correlated stripes. Synthetic calculations show that these types of signatures in the correlation coefficient matrix can result from certain combinations of cross-frequency correlated random components and diffuse field. The analysis of observed data indicates that the secondary microseismic peak around 0.15 Hz is correlated with its neighbouring frequencies, while the primary peak around 0.06 Hz is more diffuse. This suggests that the primary and secondary peaks may be associated with somewhat different physical origins. In addition, significant correlation of frequencies below that of the primary microseismic peak suggests that the very low frequencies noise is less scattered during propagation. The power spectra recorded by a station close to the edge of the Los Angeles basin is higher compared to data recorded by stations outside the basin perhaps because of enhanced basin reverberations and/or closer proximity to the ocean. This and other regional variations should be tested further using data from many more stations.

Description: An improved description of scattering and inverse scattering processes in reflection seismology may be obtained on the basis of a scattering series solution to the Helmoltz equation, which allows one to separately model primary and multiple reflections. However, the popular scattering series of Born is of limited seismic modelling value, since it is only guaranteed to converge if the global contrast is relatively small. For frequency-domain waveform modelling of realistic contrasts, some kind of renormalization may be required. The concept of renormalization is normally associated with quantum field theory, where it is absolutely essential for the treatment of infinities in connection with observable quantities. However, the renormalization program is also highly relevant for classical systems, especially when there are interaction effects that act across different length scales. In the scattering series of De Wolf, a renormalization of the Green's functions is achieved by a split of the scattering potential operator into fore- and backscattering parts; which leads to an effective reorganization and partially re-summation of the different terms in the Born series, so that their order better reflects the physics of reflection seismology. It has been demonstrated that the leading (single return) term in the De Wolf series (DWS) gives much more accurate results than the corresponding Born approximation, especially for models with high contrasts that lead to a large accumulation of phase changes in the forward direction. However, the higher order terms in the DWS that are associated with internal multiples have not been studied numerically before. In this paper, we report from a systematic numerical investigation of the convergence properties of the DWS which is based on two new operator representations of the DWS. The first operator representation is relatively similar to the original scattering potential formulation, but more global and explicit in nature. The second representation is based on the T -operator formalism from quantum scattering theory, that offers a different perspective on the interaction between up- and downgoing waves, as well as significant computational advantages (e.g. domain decomposition and fast recursive methods for one-way propagators). Our numerical results demonstrate the convergence properties of the DWS are indeed superior to those of the Born series.

Description: In this study, we develop an efficient boundary integral equation method for estimation of seismic motion in a graded medium with multiple cavities under antiplane strain conditions. This inhomogeneous and heterogeneous medium is subjected to either time-harmonic incident shear seismic waves or to body waves radiating from a point seismic source. Three different types of soil material gradient are considered: (i) density and shear modulus vary proportionally as quadratic functions of depth, but the wave velocity remains constant; (ii) the soil material is viscoelastic, with a shear modulus and density that vary with respect to the spatial coordinates in an arbitrary fashion, so that the wave velocity is both frequency and position-dependent and (iii) the soil material has position-dependent shear modulus and constant density, yielding a linear profile for the wave velocity. Three different, frequency-dependent boundary integral equation schemes are respectively developed for the aforementioned three types of graded soil materials based on: (i) Green's function for the quadratically graded elastic half-plane; (ii) a fundamental solution for the viscoelastic full-plane with position-dependent wave speed profiles and (iii) a fundamental solution for an elastic full-plane with a linearly varying wave speed profile. Next, a number of cases involving geological media with position-dependent material properties and any number of cavities of various shapes and geometry are solved in the frequency domain. The numerical results reveal the dependency of the wave fields and zones of stress concentration on the following key factors: (i) type and properties of the soil material gradient; (ii) type and characteristics of the applied seismic load; (iii) shape, position and number of cavities and (iv) interaction phenomena between the cavities and the free surface.

Description: Seismic moment tensor is one of the most important source parameters defining the earthquake dimension and style of the activated fault. Geoscientists ordinarily use moment tensor catalogues, however, few attempts have been done to assess possible impacts of moment magnitude uncertainties upon their analysis. The 2012 May 20 Emilia main shock is a representative event since it is defined in literature with a moment magnitude value ( M w ) spanning between 5.63 and 6.12. A variability of ~0.5 units in magnitude leads to a controversial knowledge of the real size of the event and reveals how the solutions could be poorly constrained. In this work, we investigate the stability of the moment tensor solution for this earthquake, studying the effect of five different 1-D velocity models, the number and the distribution of the stations used in the inversion procedure. We also introduce a 3-D velocity model to account for structural heterogeneity. We finally estimate the uncertainties associated to the computed focal planes and the obtained M w . We conclude that our reliable source solutions provide a moment magnitude that ranges from 5.87, 1-D model, to 5.96, 3-D model, reducing the variability of the literature to ~0.1. We endorse that the estimate of seismic moment from moment tensor solutions, as well as the estimate of the other kinematic source parameters, requires coming out with disclosed assumptions and explicit processing workflows. Finally and, probably more important, when moment tensor solution is used for secondary analyses it has to be combined with the same main boundary conditions (e.g. wave-velocity propagation model) to avoid conflicting results.

Description: Trade-offs between velocity and anisotropy heterogeneity complicate the interpretation of differential traveltime data and have the potential to bias isotropic tomographic models. By constructing a simple parametrisation to describe an elastic tensor with hexagonal symmetry, we find analytic solutions to the Christoffel equations in terms of fast and slow horizontal velocities that allow us to simultaneously invert differential traveltime data and splitting data from teleseismic S arrivals to recover 3-D velocity and anisotropy structure. This technique provides a constraint on the depth-extent of shallow anisotropy, otherwise absent from interpretations based on SKS splitting alone. This approach is well suited to the young Woodlark Rift, where previous studies have found strong velocity variation and substantial SKS splitting in a continental rift with relatively simple geometry. This study images a low-velocity rift axis with ≤4 per cent spreading-parallel anisotropy at 50–100 km depth that separates regions of pre-existing lithospheric fabric, indicating the synchronous development of extensional crystallographic preferred orientation and lithospheric thinning. A high-velocity slab fragment north of the rift axis is associated with strike-parallel anisotropic fast axes, similar to that seen in the shallow mantle of some subduction zones. In addition to the insights provided by the anisotropy structure, the improvement in fit to the differential traveltime data demonstrates the merit to a joint inversion that accounts for anisotropy.

Description: Short-period small magnitude seismograms mainly comprise scattered waves in the form of coda waves (the tail part of the seismogram, starting after S waves and ending when the noise prevails), spanning more than 70 per cent of the whole seismogram duration. Corresponding coda envelopes provide important information about the earth inhomogeneity, which can be stochastically modeled in terms of distribution of scatterers in a random medium. In suitable experimental conditions (i.e. high earth heterogeneity), either the two parameters describing heterogeneity (scattering coefficient), intrinsic energy dissipation (coefficient of intrinsic attenuation) or a combination of them (extinction length and seismic albedo) can be used to image Earth structures. Once a set of such parameter couples has been measured in a given area and for a number of sources and receivers, imaging their space distribution with standard methods is straightforward. However, as for finite-frequency and full-waveform tomography, the essential problem for a correct imaging is the determination of the weighting function describing the spatial sensitivity of observable data to scattering and absorption anomalies. Due to the nature of coda waves, the measured parameter couple can be seen as a weighted space average of the real parameters characterizing the rock volumes illuminated by the scattered waves. This paper uses the Monte Carlo numerical solution of the Energy Transport Equation to find approximate but realistic 2-D space-weighting functions for coda waves. Separate images for scattering and absorption based on these sensitivity functions are then compared with those obtained with commonly used sensitivity functions in an application to data from an active seismic experiment carried out at Deception Island (Antarctica). Results show that these novel functions are based on a reliable and physically grounded method to image magnitude and shape of scattering and absorption anomalies. Their extension to 3-D holds promise to improve our ability to model volcanic structures using coda waves.

Description: We investigate the relationship between seismic moment M 0 and source duration t w of microearthquakes by using high-quality seismic data recorded with a vertical borehole array installed in central Taiwan. We apply a waveform cross-correlation method to the three-component records and identify several event clusters with high waveform similarity, with event magnitudes ranging from 0.3 to 2.0. Three clusters—Clusters A, B and C—contain 11, 8 and 6 events with similar waveforms, respectively. To determine how M 0 scales with t w , we remove path effects by using a path-averaged Q . The results indicate a nearly constant t w for events within each cluster, regardless of M 0 , with mean values of t w being 0.058, 0.056 and 0.034 s for Clusters A, B and C, respectively. Constant t w , independent of M 0 , violates the commonly used scaling relation ${t_w} \propto M_0^{1/3}$ . This constant duration may arise either because all events in a cluster are hosted on the same isolated seismogenic patch, or because the events are driven by external factors of constant duration, such as fluid injections into the fault zone. It may also be related to the earthquake nucleation size.

Description: Earthquake absolute location errors which can be encountered in an underground reservoir are investigated. In such an exploitation context, earthquake hypocentre errors can have an impact on the field development and economic consequences. The approach using the state-of-the-art techniques covers both the location uncertainty and the location inaccuracy—or bias—problematics. It consists, first, in creating a 3-D synthetic seismic cloud of events in the reservoir and calculating the seismic traveltimes to a monitoring network assuming certain propagation conditions. In a second phase, the earthquakes are relocated with assumptions different from the initial conditions. Finally, the initial and relocated hypocentres are compared. As a result, location errors driven by the seismic onset time picking uncertainties and inaccuracies are quantified in 3-D. Effects induced by erroneous assumptions associated with the velocity model are also modelled. In particular, 1-D velocity model uncertainties, a local 3-D perturbation of the velocity and a 3-D geostructural model are considered. The present approach is applied to the site of Rittershoffen (Alsace, France), which is one of the deep geothermal fields existing in the Upper Rhine Graben. This example allows setting realistic scenarios based on the knowledge of the site. In that case, the zone of interest, monitored by an existing seismic network, ranges between 1 and 5 km depth in a radius of 2 km around a geothermal well. Well log data provided a reference 1-D velocity model used for the synthetic earthquake relocation. The 3-D analysis highlights the role played by the seismic network coverage and the velocity model in the amplitude and orientation of the location uncertainties and inaccuracies at subsurface levels. The location errors are neither isotropic nor aleatoric in the zone of interest. This suggests that although location inaccuracies may be smaller than location uncertainties, both quantities can have a cumulative effect. Besides, small velocity uncertainties applied to the whole 1-D profile can lead to large increase of the location uncertainties. However, local variations of the velocity field around the well may have negligible effects that would make such a feature undetectable with an absolute location method. Although the reference 1-D velocity model was built from well log data, the results show that it is not a good representative of a more realistic 3-D model including a fault and its associated block shift. The amplitude and distribution of the induced location inaccuracies are such that the positioning and the orientation of features delineated by seismicity are distorted and may be difficult to correctly interpret.

Description: Compared with a single GPS system, GPS/GLONASS observations can improve the satellite visibility, optimize the spatial geometry and improve the precise positioning performance. Although the advantage over GPS-only methods in terms of positioning is clear, the potential contributions of GPS/GLONASS to co-seismic displacement determination and the subsequent seismic source inversion still require extensive study and validation. In this paper, we first extended a temporal point positioning model from GPS-only to GPS/GLONASS observations. Using this new model, the performance of the GPS/GLONASS method for obtaining co-seismic displacements was then validated via eight outdoor experiments on a shaking table. Our result reveals that the GPS/GLONASS method provides more accurate and robust co-seismic displacements than the GPS-only observations in a non-optimal observation environment. Furthermore, as a case study, observation data recorded during the September 2015 M w 8.3 Illapel earthquake in Chile were re-processed. At some stations, obvious biases were found between the co-seismic displacements derived from GPS-only and GPS/GLONASS observations. The subsequent slip distribution inversion on a curved fault confirms that the differences in the co-seismic displacements causes differences in the inversion results and that the slip distributions of the Illapel earthquake inferred from the GPS/GLONASS observations tend to be shallower and larger.

Description: Incorporation of attenuation into the normal mode sum representations of seismic signals is commonly effected by applying perturbation theory. This is fine for weak attenuation, but problematic for stronger attenuation. In this work, modes of the anelastic medium are represented as complex superpositions of elastic eigenfunctions. For the P – SV system, a generalized eigenvalue equation for the complex eigenwavenumbers and complex coefficients used to construct the anelastic eigenfunctions is derived. The generalized eigenvalue problem for the P – SV problem is exactly linear in the eigenwavenumber at the expense of doubling the dimension. The SH problem is exactly linear in the square of the eigenwavenumber. This is in contrast to a similar standing wave problem for the earth free oscillations. Attenuation is commonly incorporated into synthetic seismogram calculations by introduction of complex frequency-dependent elastic moduli. The moduli depend nonlinearly on the frequency. The independent variable in the standing wave free oscillation problem is the frequency, which makes the eigenvalue problem nonlinear. The choice of the wavenumber as the independent variable for the travelling wave problem leads to a linear problem. The Earth model may be transversely isotropic. Compressional waves and both polarizations of shear waves ( SV, SH ) are treated.

Description: We present a new technique for deriving detailed information on seismic velocities of the subsurface material from continuous ambient noise recorded by spatially dense seismic arrays. This method uses iterative double beamforming between various subarrays to extract surface wave contributions from the ambient-noise data in complex environments with unfavourable noise-source distributions. The iterative double beamforming extraction makes it possible to retrieve large amounts of Rayleigh wave traveltime information in a wide frequency band. The method is applied to data recorded by a highly dense Nodal array with 1108 vertical geophones, centred on the damage zone of the Clark branch of the San Jacinto Fault Zone south of Anza, California. The array covers a region of ~650 x 700 m 2 , with instrument spacing of 10–30 m, and continuous recording at 500 samples s –1 over 30 d in 2014. Using this iterative double beamforming on subarrays of 25 sensors and cross-correlations between all of the station pairs, we separate surface waves from body waves that are abundant in the raw cross-correlation data. Focusing solely on surface waves, maps of traveltimes are obtained at different frequencies with unprecedented accuracy at each point of a 15-m-spacing grid. Group velocity inversions at 2–4 Hz reveal depth and lateral variations in the structural properties within and around the San Jacinto Fault Zone in the study area. This method can be used over wider frequency ranges and can be combined with other imaging techniques, such as eikonal tomography, to provide unprecedented detailed structural images of the subsurface material.

Description: Knowledge of the quality factor of near-surface materials is of fundamental interest in various applications. Attenuation can be very strong close to the surface and thus needs to be properly assessed. In recent years, several researchers have studied the retrieval of attenuation coefficients from the cross correlation of ambient seismic noise. Yet, the determination of exact amplitude information from noise-correlation functions is, in contrast to the extraction of traveltimes, not trivial. Most of the studies estimated attenuation coefficients on the regional scale and within the microseism band. In this paper, we investigate the possibility to derive attenuation coefficients from seismic noise at much shallower depths and higher frequencies (〉1 Hz). The Euroseistest area in northern Greece offers ideal conditions to study quality factor retrieval from ambient noise for different rock types. Correlations are computed between the stations of a small scale array experiment (station spacings 〈2 km) that was carried out in the Euroseistest area in 2011. We employ the correlation of the coda of the correlation (C 3 ) method instead of simple cross correlations to mitigate the effect of uneven noise source distributions on the correlation amplitude. Transient removal and temporal flattening are applied instead of 1-bit normalization in order to retain relative amplitudes. The C 3 method leads to improved correlation results (higher signal-to-noise ratio and improved time symmetry) compared to simple cross correlations. The C 3 functions are rotated from the ZNE to the ZRT system and we focus on Love wave arrivals on the transverse component and on Love wave quality factors Q L . The analysis is performed for selected stations being either situated on soft soil or on weathered rock. Phase slowness is extracted using a slant-stack method. Attenuation parameters are inferred by inspecting the relative amplitude decay of Love waves with increasing interstation distance. We observe that the attenuation coefficient and Q L can be reliably extracted for stations situated on soft soil whereas the derivation of attenuation parameters is more problematic for stations that are located on weathered rock. The results are in acceptable conformance with theoretical Love wave attenuation curves that were computed using 1-D shear wave velocity and quality factor profiles from the Euroseistest area.

Description: Large earthquakes from the intermediate-depth Vrancea seismic zone are known to produce in Bucharest ground motion characterized by predominant long periods. This phenomenon has been interpreted as the combined effect of both seismic source properties and site response of the large sedimentary basin. The thickness of the unconsolidated Quaternary deposits beneath the city is more than 200 m, the total depth of sediments is more than 1000 m. Complex basin geometry and the low seismic wave velocities of the sediments are primarily responsible for the large amplification and long duration experienced during earthquakes. For a better understanding of the geological structure under Bucharest, a number of investigations using non-invasive methods have been carried out. With the goal to analyse and extract the polarization and dispersion characteristics of the surface waves, ambient vibrations and low-magnitude earthquakes have been investigated using single station and array techniques. Love and Rayleigh dispersion curves (including higher modes), Rayleigh waves ellipticity and SH -wave fundamental frequency of resonance ( f 0 SH ) have been inverted simultaneously to estimate the shear wave velocity structure under Bucharest down to a depth of about 8 km. Information from existing borehole logs was used as prior to reduce the non-uniqueness of the inversion and to constrain the shallow part of the velocity model (〈300 m). In this study, we use data from a 35-km diameter array (the URS experiment) installed by the National Institute for Earth Physics and by the Karlsruhe Institute of Technology during 10 months in the period 2003–2004. The array consisted of 32 three-component seismological stations, deployed in the urban area of Bucharest and adjacent zones. The large size of the array and the broad-band nature of the available sensors gave us the possibility to characterize the surface wave dispersion at very low frequencies (0.05–1 Hz) using frequency–wavenumber techniques. This is essential to explore and resolve the deeper portions of the basin. The horizontal to vertical spectral ratio (H/V) curves provide important additional information about the structure and are here characterized by two major peaks. The first is attributed to the fundamental frequency of the basin, while the second can be interpreted as a mixture of the second higher mode of Rayleigh waves and other types of waves such as SH waves. This hypothesis has been verified by comparing the H/V curves with the SH -wave transfer function from the retrieved velocity structure. We could also approximate the SH transfer function with H/V ratios of earthquake recordings, providing additional verification of the robustness of the proposed velocity model. The Cretaceous bedrock depth was then inverted at each URS station from the fundamental frequency of resonance and using this model. A 3-D geophysical model for Bucharest has been constructed based on the integration of the inverted velocity profiles and the available geological information using a geographic information system.

Description: Some modern seismicity in the magnitude range of M 4 and M 6 in California and eastern North America preferentially occurs at the edges of past large ruptures. Once a large earthquake rupture has occurred, stress is concentrated at the edges of the rupture, and apparently this stress concentration can trigger earthquakes at or near the rupture edges many decades or even longer after a main shock. Furthermore, the modern M ≥ 4 earthquakes in the vicinity of a past main shock usually have the same focal mechanism as the earlier main shock. There are a number of examples of this in California and Nevada, where there is a statistically significant correlation of the locations of M ≥ 4 earthquakes and the edges of 19th and 20th century fault ruptures in M w ≥ 6.5 earthquakes. In contrast, the M ≥ 4 earthquakes near the epicentres of future ruptures in California are randomly scattered around the fault with no concentration near the ends of the future fault rupture. The concentration of earthquakes near the ends of earlier large ruptures in California becomes progressively less pronounced as the smallest magnitude in the data set is reduced from M 4.0 to M 3.0. These observations also appear to be true for intraplate regions where aftershock sequences can last millennia. The identification of modern rupture-edge M ≥ 4 aftershocks can be used to help discover where and when past strong earthquakes took place, even if there is no historical record of the main shock. This is of great importance for seismic hazard studies.

Description: Surface wave methods provide a cost effective means of developing shear wave velocity ( Vs ) profiles for applications such as dynamic site characterization and seismic site response analyses. However, the inverse problem involved in obtaining a realistic layered earth model from surface wave dispersion data is inherently ill-posed, non-linear and mix-determined, without a unique solution. When available, a priori information such as geotechnical boreholes or geologic well logs should be used to aid in constraining site-specific inversion parameters. Unfortunately, a priori information is often unavailable, particularly at significant depths, and a ‘blind analysis’ must be performed. In these situations, the analyst must decide on an appropriate number of layers and ranges for their corresponding inversion parameters (i.e. trial number of layers and ranges in their respective thicknesses, shear wave velocities, compression wave velocities and mass densities). Selection of these parameters has been shown to significantly impact the results of an inversion. This paper presents a method for conducting multiple inversions utilizing systematically varied inversion layering parametrizations in order to identify and encompass the most reasonable layered earth models for a site. Each parametrization is defined by a unique layering ratio, which represents a multiplier that systemically increases the potential thickness of each layer in the inversion parametrization based on the potential thickness of the layer directly above it. The layering ratio method is demonstrated at two sites associated with the InterPacific Project, wherein it is shown to significantly aid in selecting reasonable Vs profiles that are close representations of the subsurface. While the goal of the layering ratio inversion methodology is not necessarily to find the ‘optimal’ or ‘best’ Vs profile for a site, it may be successful at doing so for certain sites/datasets. However, the primary reason for using the layering ratio method is to find Vs profiles that realistically represent the uncertainty in Vs resulting from surface wave inversion, and to avoid selection of Vs profiles that are unrealistic and adversely influenced by the choice of inversion parametrization.

Description: The appraisal of tomographic models, of fundamental importance towards better understanding the Earth's interior, consists in analysing their resolution and covariance. The discrete theory of Backus–Gilbert, solving all at once the linear problems of model estimation and appraisal, aims at evaluating weighted averages of the true model parameters. Contrary to damped least-squares techniques, one key advantage of Backus–Gilbert inversion is that no subjective regularization is needed to remove the non-uniqueness of the model solution. Indeed, it is often possible to identify unique linear combinations of the parameters even when the parameters themselves are not uniquely defined. In other words, the non-uniqueness can be broken by averaging rather than regularizing. Over the past few decades, many authors have considered that, in addition to a high computational cost, it could be a clumsy affair in the presence of data errors to practically implement the Backus–Gilbert approach to large-scale tomographic applications. In this study, we introduce and adapt to seismic tomography the Subtractive Optimally Localized Averages (SOLA) method, an alternative Backus–Gilbert formulation which retains all its advantages, but is more computationally efficient and versatile in the explicit construction of averaging kernels. As a leitmotiv, we focus on global-scale S -wave tomography and show that the SOLA method can successfully be applied to large-scale, linear and discrete tomographic problems.

Description: The envelope broadening and the peak delay of the S -wavelet of a small earthquake with increasing travel distance are results of scattering by random velocity inhomogeneities in the earth medium. As a simple mathematical model, Sato proposed a new stochastic synthesis of the scalar wavelet envelope in 3-D von Kármán type random media when the centre wavenumber of the wavelet is in the power-law spectral range of the random velocity fluctuation. The essential idea is to split the random medium spectrum into two components using the centre wavenumber as a reference: the long-scale (low-wavenumber spectral) component produces the peak delay and the envelope broadening by multiple scattering around the forward direction; the short-scale (high-wavenumber spectral) component attenuates wave amplitude by wide angle scattering. The former is calculated by the Markov approximation based on the parabolic approximation and the latter is calculated by the Born approximation. Here, we extend the theory for the envelope synthesis of a wavelet in 2-D random media, which makes it easy to compare with finite difference (FD) simulation results. The synthetic wavelet envelope is analytically written by using the random medium parameters in the angular frequency domain. For the case that the power spectral density function of the random velocity fluctuation has a steep roll-off at large wavenumbers, the envelope broadening is small and frequency independent, and scattering attenuation is weak. For the case of a small roll-off, however, the envelope broadening is large and increases with frequency, and the scattering attenuation is strong and increases with frequency. As a preliminary study, we compare synthetic wavelet envelopes with the average of FD simulation wavelet envelopes in 50 synthesized random media, which are characterized by the RMS fractional velocity fluctuation = 0.05, correlation scale a  = 5 km and the background wave velocity V 0 = 4 km s –1 . We use the radiation of a 2 Hz Ricker wavelet from a point source. For all the cases of von Kármán order  = 0.1, 0.5 and 1, we find the synthetic wavelet envelopes are a good match to the characteristics of FD simulation wavelet envelopes in a time window starting from the onset through the maximum peak to the time when the amplitude decreases to half the peak amplitude.

Description: We present a thermodynamically based formulation for modelling dynamic rupture processes in the brittle crust using a continuum damage-breakage rheology. The model combines aspects of a continuum viscoelastic damage framework for brittle solids with a continuum breakage mechanics for granular flow within dynamically generated slip zones. The formulation accounts for the density of distributed cracking and other internal flaws in damaged rocks with a scalar damage parameter, and addresses the grain size distribution of a granular phase in the slip zone with a breakage parameter. A dynamic brittle instability is associated with a critical level of damage in the solid, leading to loss of convexity of the solid strain energy, localization and transition to a granular phase associated with lower energy level. The continuum damage-breakage rheology model treats the localization to a slip zone at the onset of dynamic rupture and post-failure recovery process as phase transitions between solid and granular states. The model generates sub- and supershear rupture velocities and pulse-type ruptures seen also in frictional models, and additional important features such as strong dynamic changes of volumetric strain near the rupture front and diversity of nucleation mechanisms. The propagation of rupture front and slip accumulation at a point are correlated with sharp dynamic dilation followed by a gradual decay to a level associated with the final volumetric change associated with the granular phase transition in the slipping zone. The local brittle failure process associated with the solid–granular transition is expected to produce isotropic radiation in addition to the deviatoric terms. The framework significantly extends the ability to model brittle processes in complex geometrical structures and allows analysing the roles of gouge thickness and other parameters on nucleation, rupture and radiation characteristics.

Description: Numerical modeling of seismic waves in heterogeneous porous reservoir rocks is an important tool for the interpretation of seismic surveys in reservoir engineering. We apply globally optimal implicit staggered-grid finite differences (FD) to model 2-D wave propagation in heterogeneous poroelastic media at a low-frequency range (〈10 kHz). We validate the numerical solution by comparing it to an analytical-transient solution obtaining clear seismic wavefields including fast P and slow P and S waves (for a porous media saturated with fluid). The numerical dispersion and stability conditions are derived using von Neumann analysis, showing that over a wide range of porous materials the Courant condition governs the stability and this optimal implicit scheme improves the stability of explicit schemes. High-order explicit FD can be replaced by some lower order optimal implicit FD so computational cost will not be as expensive while maintaining the accuracy. Here, we compute weights for the optimal implicit FD scheme to attain an accuracy of = 10 –8 . The implicit spatial differentiation involves solving tridiagonal linear systems of equations through Thomas’ algorithm.

Description: S -wave receiver functions (SRFs) are used to investigate crustal and upper-mantle structure beneath several ice-covered areas of Antarctica. Moho S-to-P (Sp) arrivals are observed at ~6–8 s in SRF stacks for stations in the Gamburtsev Mountains (GAM) and Vostok Highlands (VHIG), ~5–6 s for stations in the Transantarctic Mountains (TAM) and the Wilkes Basin (WILK), and ~3–4 s for stations in the West Antarctic Rift System (WARS) and the Marie Byrd Land Dome (MBLD). A grid search is used to model the Moho Sp conversion time with Rayleigh wave phase velocities from 18 to 30 s period to estimate crustal thickness and mean crustal shear wave velocity. The Moho depths obtained are between 43 and 58 km for GAM, 36 and 47 km for VHIG, 39 and 46 km for WILK, 39 and 45 km for TAM, 19 and 29 km for WARS and 20 and 35 km for MBLD. SRF stacks for GAM, VHIG, WILK and TAM show little evidence of Sp arrivals coming from upper-mantle depths. SRF stacks for WARS and MBLD show Sp energy arriving from upper-mantle depths but arrival amplitudes do not rise above bootstrapped uncertainty bounds. The age and thickness of the crust is used as a heat flow proxy through comparison with other similar terrains where heat flow has been measured. Crustal structure in GAM, VHIG and WILK is similar to Precambrian terrains in other continents where heat flow ranges from ~41 to 58 mW m –2 , suggesting that heat flow across those areas of East Antarctica is not elevated. For the WARS, we use the Cretaceous Newfoundland–Iberia rifted margins and the Mesozoic-Tertiary North Sea rift as tectonic analogues. The low-to-moderate heat flow reported for the Newfoundland–Iberia margins (40–65 mW m –2 ) and North Sea rift (60–85 mW m –2 ) suggest that heat flow across the WARS also may not be elevated. However, the possibility of high heat flow associated with localized Cenozoic extension or Cenozoic-recent magmatic activity in some parts of the WARS cannot be ruled out.

Description: The 2014 April 1, M w 8.1 Iquique (Chile) earthquake struck in the Northern Chile seismic gap. With a rupture length of less than 200 km, it left unbroken large segments of the former gap. Early studies were able to model the main rupture features but results are ambiguous with respect to the role of aseismic slip and left open questions on the remaining hazard at the Northern Chile gap. A striking observation of the 2014 earthquake has been its extensive preparation phase, with more than 1300 events with magnitude above M L 3, occurring during the 15 months preceding the main shock. Increasing seismicity rates and observed peak magnitudes accompanied the last three weeks before the main shock. Thanks to the large data sets of regional recordings, we assess the precursor activity, compare foreshocks and aftershocks and model rupture preparation and rupture effects. To tackle inversion challenges for moderate events with an asymmetric network geometry, we use full waveforms techniques to locate events, map the seismicity rate and derive source parameters, obtaining moment tensors for more than 300 events (magnitudes M w 4.0–8.1) in the period 2013 January 1–2014 April 30. This unique data set of fore- and aftershocks is investigated to distinguish rupture process models and models of strain and stress rotation during an earthquake. Results indicate that the spatial distributions of foreshocks delineated the shallower part of the rupture areas of the main shock and its largest aftershock, well matching the spatial extension of the aftershocks cloud. Most moment tensors correspond to almost pure double couple thrust mechanisms, consistent with the slab orientation. Whereas no significant differences are observed among thrust mechanisms in different areas, nor among thrust foreshocks and aftershocks, the early aftershock sequence is characterized by the presence of normal fault mechanisms, striking parallel to the trench but dipping westward. These events likely occurred in the shallow wedge structure close to the slab interface and are consequence of the increased extensional stress in this region after the largest events. The overall stress inversion result suggests a minor stress rotation after the main shock, but a significant release of the deviatoric stress. The temporal change in the distribution of focal mechanisms can also be explained in terms of the spatial heterogeneity of the stress field: under such interpretation, the potential of a large megathrust earthquake breaking a larger segment offshore Northern Chile remains high.

Description: We analyse the performance of a higher order accurate staggered viscoelastic time-domain finite-difference method, in which the staggered Adams–Bashforth (ABS) third-order and fourth-order accurate time integrators are used for temporal discretization. ABS is a multistep method that uses previously calculated wavefields to increase the order of accuracy in time. The analysis shows that the numerical dispersion is much lower than that of the widely used second-order leapfrog method. Numerical dissipation is introduced by the ABS method which is significantly smaller for fourth-order than third-order accuracy. In 1-D and 3-D simulation experiments, we verify the convincing improvements of simulation accuracy of the fourth-order ABS method. In a realistic elastic 3-D scenario, the computing time reduces by a factor of approximately 2.4, whereas the memory requirements increase by approximately a factor of 2.2. The ABS method thus provides an alternative strategy to increase the simulation accuracy in time by investing computer memory instead of computing time.

Description: Synthetic dispersion curves are here computed in the frame of an ambient-vibration full-wavefield model, which relies on the description of both ambient-vibration ground displacement and its sources as stochastic fields defined on the Earth's surface, stationary in time and homogeneous in space. In this model, previously developed for computing synthetic Horizontal-to-Vertical Spectral Ratio curves, the power spectral density function and the spatial autocorrelation of the displacement are naturally described as functions of the power spectral density function of the generating forces and of the subsoil properties ( via the relevant Green's function), by also accounting for spatial correlation of these forces. Dispersion curves are computed from the displacement power spectral density function and from the spatial autocorrelation according with the well-known f-k and SPAC techniques, respectively. Two examples illustrate the way this new ambient-vibration model works, showing its possible use in better understanding the role of the surface waves in forming the dispersion curves, as well as its capability to capture some remarkable experimental findings.

Description: The Abruzzi region is located in the Central Apennines Neogene fold-and-thrust belt and has one of the highest seismogenic potential in Italy, with high and diffuse crustal seismicity related to NE–SW oriented extension. In this study, we investigate the detailed spatial variation in shear wave splitting providing high-resolution anisotropic structure beneath the L'Aquila region. To accomplish this, we performed a systematic analysis of crustal anisotropic parameters: fast polarization direction ( ) and delay time ( t ). We benefit from the dense coverage of seismic stations operating in the area and from a catalogue of several accurate earthquake locations of the 2009 L'Aquila seismic sequence, related to the M w 6.1 2009 L'Aquila main shock, to describe in detail the geometry of the anisotropic volume around the active faults that ruptured. The spatial variations both in and t suggest a complex anisotropic structure beneath the region caused by a combination of both structural- and stress-induced mechanisms. The average is NNW–SSE oriented (N141°), showing clear similarity both with the local fault strike and the SH max . In the central part of the study area fast axes are oriented NW–SE, while moving towards the northeastern and northwestern sectors the fast directions clearly diverge from the general trend of NW–SE and rotate accordingly to the local fault strikes. The above-mentioned fault-parallel distribution suggests that the observed anisotropy is mostly controlled by the local fault-related structure. Toward the southeast fast directions become orthogonal both to strike of the local mapped faults and to the SH max . Here, are predominantly oriented NE–SW; we interpret this orientation as due to the presence of a highly fractured and overpressurized rock volume which should be responsible of the 90° flips in and the increase in t . Another possible mechanism for NE–SW orientation of in the southeastern sector could be ascribed to the presence of a buried, deep NE–SW oriented fault system. t , both unnormalized and normalized, does not show any clear evidence of increasing with increasing depth indicating that the anisotropy is confined primarily to the shallower crustal layers (~10 km depth). Interpolating t show that higher values are found at the edges of the main patches of the rupture related to the 2009 main shock, while lower values are limited in the central part of the fault plane, where the coseismic slip was higher. We infer that in the areas surrounding the ruptured region, lateral variations in material properties caused overpressurized fluid conditions, while within the main shock ruptured area, high energy released produced an open crack system such that overpressurization was not possible.

Description: This paper presents a set of deterministic 3-D ground motion simulations for the greater metropolitan area of Medellín in the Aburrá Valley, an earthquake-prone region of the Colombian Andes that exhibits moderate-to-strong topographic irregularities. We created the velocity model of the Aburrá Valley region (version 1) using the geological structures as a basis for determining the shear wave velocity. The irregular surficial topography is considered by means of a fictitious domain strategy. The simulations cover a 50 x 50 x 25 km 3 volume, and four M w = 5 rupture scenarios along a segment of the Romeral fault, a significant source of seismic activity in Colombia. In order to examine the sensitivity of ground motion to the irregular topography and the 3-D effects of the valley, each earthquake scenario was simulated with three different models: (i) realistic 3-D velocity structure plus realistic topography, (ii) realistic 3-D velocity structure without topography, and (iii) homogeneous half-space with realistic topography. Our results show how surface topography affects the ground response. In particular, our findings highlight the importance of the combined interaction between source-effects, source-directivity, focusing, soft-soil conditions, and 3-D topography. We provide quantitative evidence of this interaction and show that topographic amplification factors can be as high as 500 per cent at some locations. In other areas within the valley, the topographic effects result in relative reductions, but these lie in the 0–150 per cent range.

Description: Wave-induced local fluid flow is known as a key mechanism to explain the intrinsic wave dissipation in fluid-saturated rocks. Understanding the relationship between the acoustic properties of rocks and fluid patch distributions is important to interpret the observed seismic wave phenomena. A triple-layer patchy (TLP) model is proposed to describe the P -wave dissipation process in a double-porosity media saturated with two immiscible fluids. The double-porosity rock consists of a solid matrix with unique host porosity and inclusions which contain the second type of pores. Two immiscible fluids are considered in concentric spherical patches, where the inner pocket and the outer sphere are saturated with different fluids. The kinetic and dissipation energy functions of local fluid flow (LFF) in the inner pocket are formulated through oscillations in spherical coordinates. The wave propagation equations of the TLP model are based on Biot's theory and the corresponding Lagrangian equations. The P -wave dispersion and attenuation caused by the Biot friction mechanism and the local fluid flow (related to the pore structure and the fluid distribution) are obtained by a plane-wave analysis from the Christoffel equations. Numerical examples and laboratory measurements indicate that P -wave dispersion and attenuation are significantly influenced by the spatial distributions of both, the solid heterogeneity and the fluid saturation distribution. The TLP model is in reasonably good agreement with White's and Johnson's models. However, differences in phase velocity suggest that the heterogeneities associated with double-porosity and dual-fluid distribution should be taken into account when describing the P -wave dispersion and attenuation in partially saturated rocks.

Description: Seismic ambient noise tomography is applied to central and southern Mozambique, located in the tip of the East African Rift (EAR). The deployment of MOZART seismic network, with a total of 30 broad-band stations continuously recording for 26 months, allowed us to carry out the first tomographic study of the crust under this region, which until now remained largely unexplored at this scale. From cross-correlations extracted from coherent noise we obtained Rayleigh wave group velocity dispersion curves for the period range 5–40 s. These dispersion relations were inverted to produce group velocity maps, and 1-D shear wave velocity profiles at selected points. High group velocities are observed at all periods on the eastern edge of the Kaapvaal and Zimbabwe cratons, in agreement with the findings of previous studies. Further east, a pronounced slow anomaly is observed in central and southern Mozambique, where the rifting between southern Africa and Antarctica created a passive margin in the Mesozoic, and further rifting is currently happening as a result of the southward propagation of the EAR. In this study, we also addressed the question concerning the nature of the crust (continental versus oceanic) in the Mozambique Coastal Plains (MCP), still in debate. Our data do not support previous suggestions that the MCP are floored by oceanic crust since a shallow Moho could not be detected, and we discuss an alternative explanation for its ocean-like magnetic signature. Our velocity maps suggest that the crystalline basement of the Zimbabwe craton may extend further east well into Mozambique underneath the sediment cover, contrary to what is usually assumed, while further south the Kaapval craton passes into slow rifted crust at the Lebombo monocline as expected. The sharp passage from fast crust to slow crust on the northern part of the study area coincides with the seismically active NNE-SSW Urema rift, while further south the Mazenga graben adopts an N-S direction parallel to the eastern limit of the Kaapvaal craton. We conclude that these two extensional structures herald the southward continuation of the EAR, and infer a structural control of the transition between the two types of crust on the ongoing deformation.

Description: Seismic wavefields interact with low-velocity fault damage zones. Waveforms of ballistic fault zone head waves, trapped waves, reflected waves and signatures of trapped noise can provide important information on structural and mechanical fault zone properties. Here we extend the class of observable fault zone waves and reconstruct in-fault reverberations or multiples in a strike-slip faulting environment. Manifestations of the reverberations are significant, consistent wave fronts in the coda of cross-correlation functions that are obtained from scattered earthquake waveforms and seismic noise recorded by a linear fault zone array. The physical reconstruction of Green's functions is evident from the high similarity between the signals obtained from the two different scattered wavefields. Modal partitioning of the reverberation wavefield can be tuned using different data normalization techniques. The results imply that fault zones create their own ambiance, and that the here reconstructed reverberations are a key seismic signature of wear zones. Using synthetic waveform modelling we show that reverberations can be used for the imaging of structural units by estimating the location, extend and magnitude of lateral velocity contrasts. The robust reconstruction of the reverberations from noise records suggests the possibility to resolve the response of the damage zone material to various external and internal loading mechanisms.

Description: The Lower Rhine Embayment in Central Europe hosts a rift system that has very low deformation rates. The faults in this area have slip rates of less than 0.1 mm yr –1 , which does not allow to investigate ongoing tectonic deformation with geodetic techniques, unless they cover very long time spans. Instrumental seismicity does only cover a small fraction of the very long earthquake recurrence intervals of several thousands of years. Palaeoseismological studies are needed to constrain slip rates and the earthquake history of such faults. Destructive earthquakes are rare in the study area, but did occur in historic times. In 1755/1756, a series of strong earthquakes caused significant destruction in the city of Düren (Germany) and the surrounding areas. In this study we document palaeoseismological data from the nearby Rurrand Fault. In contrast to earlier studies on the same fault, we found evidence for a surface rupturing earthquake in the Holocene, and we identified at least one more surface rupturing event. Our study shows that the Rurrand Fault currently accommodates deformation in earthquakes rather than by creeping. The coseismic offsets were determined to be between less than 0.5 m per event. We assign maximum possible magnitudes of M w 5.9–6.8 for the Rurrand Fault and a slip rate of at least 0.02–0.03 mm yr –1 for the last ~130–50 kyr. The surface ruptures did not occur at the main fault trace that has a clear morphological expression due to older tectonic motions, but on a younger fault strand in the hanging wall of the main fault. Terrain analyses based on 1 m resolution airborne LiDAR data have been used to image the subtle morphological expression of this young fault zone. Georadar and electric resistivity tomography were applied to image the fault zone at depth and to test if these shallow geophysical methods can be used to identify and trace the fault zone. Georadar failed to produce reliable results, but geoelectrics were successfully applied and allowed us to retrieve slip rate estimates. Our results indicate that the Düren 1755/1756 earthquakes did not produce surface ruptures at the Rurrand Fault, either because they did not rupture the surface at all, or because they occurred at another, neighbouring fault.

Description: We have measured radial anisotropy in the crust beneath the northern part of North China by jointly inverting Rayleigh and Love wave phase velocities at periods less than 35 s from 14 months of ambient noise data recorded by 222 broad-band seismic stations. We also estimate the azimuthal anisotropy of phase velocity from Rayleigh wave data. The fast direction of azimuthal anisotropy varies with periods, NE–SW orientation at short and intermediate periods (10–16 s) and NW–SE orientation at periods larger than 20 s. The NE–SW oriented fast direction of azimuthal anisotropy may be related to the fossilized structural fabrics due to the compression during the Indosinian orogeny from late Palaeozoic to middle Mesozoic. The NW–SE trend of anisotropic fabric in the lower crust and uppermost mantle is probably associated with the later lithospheric extension. The observed radial anisotropy also shows a two-layer feature, negative radial anisotropy (Vsh 〈 Vsv) in the upper crust and positive (Vsh 〉 Vsv) in the middle-lower crust. The compressional tectonics from late Palaeozoic to middle Mesozoic may cause crustal materials align vertically throughout the crust. This vertical fabric could make Vsh slower than Vsv. However, the lithospheric extension in the late Mesozoic to Cenozoic time could overprint the older fabric in the middle and lower crust by magma intrusion and underplating. Horizontal alignment of the material or intruded melt sills due to the extension probably produce the observed strong positive radial anisotropy in the middle and lower crust.

Description: Using active and passive seismology data we derive a shear ( S ) wave velocity model and a Poisson's ratio () model across the Chilean convergent margin along a profile at 38°15'S, where the M w 9.5 Valdivia earthquake occurred in 1960. The derived S -wave velocity model was constructed using three independently obtained velocity models that were merged together. In the upper part of the profile (0–2 km depth), controlled source data from explosions were used to obtain an S -wave traveltime tomogram. For the middle part (2–20 km depth), data from a temporary seismology array were used to carry out a dispersion analysis. The resulting dispersion curves were used to obtain a 3-D S -wave velocity model. In the lower part (20–75 km depth, depending on the longitude), an already existent local earthquake tomographic image was merged with the other two sections. This final S -wave velocity model and already existent compressional ( P ) wave velocity models along the same transect allowed us to obtain a Poisson's ratio model. The results of this study show that the velocities and Poisson's ratios in the continental crust of this part of the Chilean convergent margin are in agreement with geological features inferred from other studies and can be explained in terms of normal rock types. There is no requirement to call on the existence of measurable amounts of present-day fluids, in terms of seismic velocities, above the plate interface in the continental crust of the Coastal Cordillera and the Central Valley in this part of the Chilean convergent margin. This is in agreement with a recent model of water being transported down and released from the subduction zone.

Description: In the estimation of spatiotemporal slip models, kinematic source inversions using Akaike's Bayesian Information Criterion (ABIC) and the multiple-time-window method have often been used. However, there are cases in which conventional ABIC-based source inversions do not work well in the determination of hyperparameters when a non-negative slip constraint is used. In order to overcome this problem, a new source inversion method was developed in this study. The new method introduces a fully Bayesian method into the kinematic multiple-time-window source inversion. The multiple-time-window method is one common way of parametrizing a source time function and is highly flexible in terms of the shape of the source time function. The probability distributions of model parameters and hyperparameters can be directly obtained by using the Markov chain Monte Carlo method. These probability distributions are useful for simply evaluating the uniqueness and reliability of the derived model, which is another advantage of a fully Bayesian method. This newly developed source inversion method was applied to the 2011 Ibaraki-oki, Japan, earthquake ( M w 7.9) to demonstrate its usefulness. It was demonstrated that the problem with using the conventional ABIC-based source inversion method for hyperparameter determination appeared in the spatiotemporal source inversion of this event and that the newly developed source inversion could overcome this problem.

Description: In a set of problems ranging from 4-D seismic to salt boundary estimation, updates to the velocity model often have a highly localized nature. Numerical techniques for these applications such as full-waveform inversion (FWI) require an estimate of the wavefield to compute the model updates. When dealing with localized problems, it is wasteful to compute these updates in the global domain, when we only need them in our region of interest. This paper introduces a local solver that generates forward and adjoint wavefields which are, to machine precision, identical to those generated by a full-domain solver evaluated within the region of interest. This means that the local solver computes all interactions between model updates within the region of interest and the inhomogeneities in the background model outside. Because no approximations are made in the calculation of the forward and adjoint wavefields, the local solver can compute the identical gradient in the region of interest as would be computed by the more expensive full-domain solver. In this paper, the local solver is used to efficiently generate the FWI gradient at the boundary of a salt body. This gradient is then used in a level set method to automatically update the salt boundary.

Description: Full waveform inversion using the conventional L 2 distance to measure the misfit between seismograms is known to suffer from cycle skipping. An alternative strategy is proposed in this study, based on a measure of the misfit computed with an optimal transport distance. This measure allows to account for the lateral coherency of events within the seismograms, instead of considering each seismic trace independently, as is done generally in full waveform inversion. The computation of this optimal transport distance relies on a particular mathematical formulation allowing for the non-conservation of the total energy between seismograms. The numerical solution of the optimal transport problem is performed using proximal splitting techniques. Three synthetic case studies are investigated using this strategy: the Marmousi 2 model, the BP 2004 salt model, and the Chevron 2014 benchmark data. The results emphasize interesting properties of the optimal transport distance. The associated misfit function is less prone to cycle skipping. A workflow is designed to reconstruct accurately the salt structures in the BP 2004 model, starting from an initial model containing no information about these structures. A high-resolution P -wave velocity estimation is built from the Chevron 2014 benchmark data, following a frequency continuation strategy. This estimation explains accurately the data. Using the same workflow, full waveform inversion based on the L 2 distance converges towards a local minimum. These results yield encouraging perspectives regarding the use of the optimal transport distance for full waveform inversion: the sensitivity to the accuracy of the initial model is reduced, the reconstruction of complex salt structure is made possible, the method is robust to noise, and the interpretation of seismic data dominated by reflections is enhanced.

Description: This paper concerns the numerical modelling of time-domain mechanical waves in viscoelastic media based on a generalized Zener model. To do so, classically in the literature relaxation mechanisms are introduced, resulting in a set of the so-called memory variables and thus in large computational arrays that need to be stored. A challenge is thus to accurately mimic a given attenuation law using a minimal set of relaxation mechanisms. For this purpose, we replace the classical linear approach of Emmerich & Korn with a nonlinear optimization approach with constraints of positivity. We show that this technique is more accurate than the linear approach. Moreover, it ensures that physically meaningful relaxation times that always honour the constraint of decay of total energy with time are obtained. As a result, these relaxation times can always be used in a stable way in a modelling algorithm, even in the case of very strong attenuation for which the classical linear approach may provide some negative and thus unusable coefficients.

Description: Fractures are a source of extra compliance in the rock mass. The mechanical properties of the rock matrix as well as the propagation of seismic waves inside the rock medium are dependent on the magnitude of roughness and offset between the imperfect fracture interfaces. Fracture compliance can estimate the degree of contact between fracture faces, type of fluid filling the fracture and the fracture roughness. To characterize these fracture properties, compliance ratio, known by the ratio of normal-to-shear compliance, can be a potential tool in the subsurface studies to improve the well layout design. The focus of this study is to illustrate how the compliance ratio of a rough fracture, with or without the offset between the fracture faces, can diverge from the compliance ratio of a fracture with smooth interface. Quasi-static and dynamic methods are two common ways to calculate the compliance. The former calculates the compliance by measuring the change in the displacement with the applied stress, while the latter estimates the compliance through monitoring the changes in propagation of seismic waves. To compare the compliance ratios of fractures with imperfect and smooth interfaces in an infinite medium, a numerical finite-element model is built in commercial finite-element software. The imperfect interface of the fracture is modeled with saw-tooth-like structures where they can be partially or fully in contact. The defined saw-tooth-like structures of contact asperities impose an in-plane asymmetry in the shear direction. This asymmetry causes two different values for the compliance in shear direction, known as the soft and stiff shear compliance. Our numerical simulations suggest the increase in the degree of contact between the fracture faces increases the compliance ratio in the stiff direction more than the soft direction. The compliance ratio of the fracture with the imperfect interface is larger than the compliance ratio of the smooth fracture. We suggest that the interlocking and riding up effects at the fracture interface may explain our findings in this study.

Description: Ambient noise correlation method has been widely used to extract surface waves and tomography. The extraction of body waves has been very limited, but recent reports have suggested promises for deep incident waves. Here we report our first observations of triplicated PKP phases (important phases for studying the Earth's core) and confirm observations of other body-wave core phases from noise correlations. We use dense seismic arrays in South America and China Regional Seismic Networks at distances from 145° to the antipode. We can clearly observe different PKP branches (df, bc and ab) in stacks of the station–station correlations. Both ambient noise and earthquake coda contribute to PKP phases. However, the contributions vary with frequency and with body-wave phases. At shorter periods (5–20 s), three branches of PKP (df, bc and ab) can be extracted from ambient noise and the ab phase from earthquake coda. At longer periods (15–50 s), earthquake coda are effective in generating the df branch, but not the ab branch. The generation of the PKIKP phase (df branch) from earthquake coda does not depend on earthquake focal mechanisms or focal depths. However, earthquakes far from the stations contribute more than events closer by. The best coda window is around 10 000–40 000 s and the best magnitude threshold is M w greater than 6.8 or 6.9. The observation of triplicated PKP branches from noise correlations provides a new type of data for studying the Earth's deep interior, in particularly the inner core anisotropy, which overcomes some of the limitations of traditional earthquake-based studies (such as limited source distributions and source location errors).

Description: The project Seismic Hazard Harmonization in Europe (SHARE), completed in 2013, presents significant improvements over previous regional seismic hazard modeling efforts. The Global Strain Rate Map v2.1, sponsored by the Global Earthquake Model Foundation and built on a large set of self-consistent geodetic GPS velocities, was released in 2014. To check the SHARE seismic source models that were based mainly on historical earthquakes and active fault data, we first evaluate the SHARE historical earthquake catalogues and demonstrate that the earthquake magnitudes are acceptable. Then, we construct an earthquake potential model using the Global Strain Rate Map data. SHARE models provided parameters from which magnitude–frequency distributions can be specified for each of 437 seismic source zones covering most of Europe. Because we are interested in proposed magnitude limits, and the original zones had insufficient data for accurate estimates, we combine zones into five groups according to SHARE's estimates of maximum magnitude. Using the strain rates, we calculate tectonic moment rates for each group. Next, we infer seismicity rates from the tectonic moment rates and compare them with historical and SHARE seismicity rates. For two of the groups, the tectonic moment rates are higher than the seismic moment rates of the SHARE models. Consequently, the rates of large earthquakes forecast by the SHARE models are lower than those inferred from tectonic moment rate. In fact, the SHARE models forecast higher seismicity rates than the historical rates, which indicate that the authors of SHARE were aware of the potentially higher seismic activities in the zones. For one group, the tectonic moment rate is lower than the seismic moment rates forecast by the SHARE models. As a result, the rates of large earthquakes in that group forecast by the SHARE model are higher than those inferred from tectonic moment rate, but lower than what the historical data show. For the other two groups, the seismicity rates from tectonic moment rate, historical data and SHARE models are consistent. For four groups, the maximum magnitudes used by SHARE are fairly consistent with the probable maximum magnitudes inferred from tectonic strain rates. This study demonstrates that: (1) the strain-rate data are useful for constraining seismicity rates and magnitude limits; and (2) SHARE seismic source models and historical earthquake catalogues are satisfactory.

Description: We present a method for estimating ratios of P and S waves corner frequencies ( R cf ) and earthquake strain drops by joint analysis of P and S source spectra of neighbouring groups of events. The method is applied systematically to data generated by ~9000 earthquakes around the Karadere segment of the North Anatolian Fault Zone. The results indicate several regions that produce consistently R cf values higher (e.g. 〉2) than expected from classical earthquake source models. These are associated generally with fault sections having strong geometrical heterogeneities, shallow depth sections and/or locations without large pre-existing surface trace. Earthquake ruptures in such regions are likely to generate significant rock damage and tensile components of faulting. To assess whether the observed high R cf values are produced by enriched high frequency P waves, reduced high frequency S waves or both, we compare the associated P and S spectra with mean/median results. The analysis suggests that the high R cf values of shallow events (depth 〈4 km) are generated primarily by reduced high frequency S radiation, and that the contribution from elevated high frequency P radiation increases with depth and proximity to geometrical complexities. The results highlight the importance of considering carefully the existence of some volumetric source components in earthquake rupture processes.

Description: Fractures are a source of extra compliance in the rock mass. The mechanical properties of the rock matrix as well as the propagation of seismic waves inside the rock medium are dependent on the magnitude of roughness and offset between the imperfect fracture interfaces. Fracture compliance can estimate the degree of contact between fracture faces, type of fluid filling the fracture and the fracture roughness. To characterize these fracture properties, compliance ratio, known by the ratio of normal-to-shear compliance, can be a potential tool in the subsurface studies to improve the well layout design. The focus of this study is to illustrate how the compliance ratio of a rough fracture, with or without the offset between the fracture faces, can diverge from the compliance ratio of a fracture with smooth interface. Quasi-static and dynamic methods are two common ways to calculate the compliance. The former calculates the compliance by measuring the change in the displacement with the applied stress, while the latter estimates the compliance through monitoring the changes in propagation of seismic waves. To compare the compliance ratios of fractures with imperfect and smooth interfaces in an infinite medium, a numerical finite-element model is built in commercial finite-element software. The imperfect interface of the fracture is modeled with saw-tooth-like structures where they can be partially or fully in contact. The defined saw-tooth-like structures of contact asperities impose an in-plane asymmetry in the shear direction. This asymmetry causes two different values for the compliance in shear direction, known as the soft and stiff shear compliance. Our numerical simulations suggest the increase in the degree of contact between the fracture faces increases the compliance ratio in the stiff direction more than the soft direction. The compliance ratio of the fracture with the imperfect interface is larger than the compliance ratio of the smooth fracture. We suggest that the interlocking and riding up effects at the fracture interface may explain our findings in this study.

Description: Full waveform inversion using the conventional L 2 distance to measure the misfit between seismograms is known to suffer from cycle skipping. An alternative strategy is proposed in this study, based on a measure of the misfit computed with an optimal transport distance. This measure allows to account for the lateral coherency of events within the seismograms, instead of considering each seismic trace independently, as is done generally in full waveform inversion. The computation of this optimal transport distance relies on a particular mathematical formulation allowing for the non-conservation of the total energy between seismograms. The numerical solution of the optimal transport problem is performed using proximal splitting techniques. Three synthetic case studies are investigated using this strategy: the Marmousi 2 model, the BP 2004 salt model, and the Chevron 2014 benchmark data. The results emphasize interesting properties of the optimal transport distance. The associated misfit function is less prone to cycle skipping. A workflow is designed to reconstruct accurately the salt structures in the BP 2004 model, starting from an initial model containing no information about these structures. A high-resolution P -wave velocity estimation is built from the Chevron 2014 benchmark data, following a frequency continuation strategy. This estimation explains accurately the data. Using the same workflow, full waveform inversion based on the L 2 distance converges towards a local minimum. These results yield encouraging perspectives regarding the use of the optimal transport distance for full waveform inversion: the sensitivity to the accuracy of the initial model is reduced, the reconstruction of complex salt structure is made possible, the method is robust to noise, and the interpretation of seismic data dominated by reflections is enhanced.

Description: We present a novel earthquake location method using acoustic wave-equation-based traveltime inversion. The linear relationship between the location perturbation ( t 0 , x s ) and the resulting traveltime residual t of a particular seismic phase, represented by the traveltime sensitivity kernel K ( t 0 , x s ) with respect to the earthquake location ( t 0 , x s ), is theoretically derived based on the adjoint method. Traveltime sensitivity kernel K ( t 0 , x s ) is formulated as a convolution between the forward and adjoint wavefields, which are calculated by numerically solving two acoustic wave equations. The advantage of this newly derived traveltime kernel is that it not only takes into account the earthquake–receiver geometry but also accurately honours the complexity of the velocity model. The earthquake location is obtained by solving a regularized least-squares problem. In 3-D realistic applications, it is computationally expensive to conduct full wave simulations. Therefore, we propose a 2.5-D approach which assumes the forward and adjoint wave simulations within a 2-D vertical plane passing through the earthquake and receiver. Various synthetic examples show the accuracy of this acoustic wave-equation-based earthquake location method. The accuracy and efficiency of the 2.5-D approach for 3-D earthquake location are further verified by its application to the 2004 Big Bear earthquake in Southern California.

Description: Seismic interferometry applied to 120 hr of railroad traffic recorded by an array of vertical component seismographs along a railway within the Rio Grande rift has recovered surface and body waves characteristic of the geology beneath the railway. Linear and hyperbolic arrivals are retrieved that agree with surface (Rayleigh), direct and reflected P waves observed by nearby conventional seismic surveys. Train-generated Rayleigh waves span a range of frequencies significantly higher than those recovered from typical ambient noise interferometry studies. Direct P -wave arrivals have apparent velocities appropriate for the shallow geology of the survey area. Significant reflected P -wave energy is also present at relatively large offsets. A common midpoint stack produces a reflection image consistent with nearby conventional reflection data. We suggest that for sources at the free surface (e.g. trains) increasing the aperture of the array to record wide angle reflections, in addition to longer recording intervals, might allow the recovery of deeper geological structure from railroad traffic. Frequency–wavenumber analyses of these recordings indicate that the train source is symmetrical (i.e. approaching and receding) and that deeper refracted energy is present although not evident in the time-offset domain. These results confirm that train-generated vibrations represent a practical source of high-resolution subsurface information, with particular relevance to geotechnical and environmental applications.

Description: An approach to depth migration, based on an integral representation of seismic data, that is, wavefields recorded on the boundary, is presented in terms of Poincaré wavelets. Each wavelet is taken as a boundary datum for a high-frequency asymptotic solution of the wave equation. This solution, which we call the quasiphoton or the Gaussian wave packet, decreases in a Gaussian manner away from a point running along a ray that is launched from the surface. The deformation of the propagating packet is taken into account in the migration algorithm. A numerical example of zero-offset migration with synthetic seismograms calculated for the 2-D SEG/EAGE salt model is presented. The result, which uses only 3.9 per cent of the total number of coefficients, is a satisfactory image, with a threshold of 0.75 per cent.

Description: The motion of a self-gravitating hyperelastic body is described through a time-dependent mapping from a reference body into physical space, and its material properties are determined by a referential density and strain-energy function defined relative to the reference body. Points within the reference body do not have a direct physical meaning, but instead act as particle labels that could be assigned in different ways. We use Hamilton's principle to determine how the referential density and strain-energy functions transform when the particle labels are changed, and describe an associated ‘particle relabelling symmetry’. We apply these results to linearized elastic wave propagation and discuss their implications for seismological inverse problems. In particular, we show that the effects of boundary topography on elastic wave propagation can be mapped exactly into volumetric heterogeneity while preserving the form of the equations of motion. Several numerical calculations are presented to illustrate our results.

Description: Scattering due to small-scale heterogeneities in the rock formation surrounding a wellbore can significantly change the acoustic waveform from a logging measurement which in turn can be used to characterize the formation heterogeneities. This study simulates the elastic heterogeneity scattering in monopole and dipole acoustic logging and analyse the resulting effects on the waveforms. The results show that significant coda waves are generated in both monopole and dipole waveforms and the dipole coda is dominated by S -to- S scattering, which can be effectively utilized to diagnose the heterogeneity in the rock formation. The coda wave modelling and analysis were used to characterize dipole acoustic data logged before and after fracturing a reservoir interval, with significant coda wave in the after-fracturing data indicating fracturing-induced heterogeneous property change in the rock volume surrounding the borehole.

Description: We present a novel earthquake location method using acoustic wave-equation-based traveltime inversion. The linear relationship between the location perturbation ( t 0 , x s ) and the resulting traveltime residual t of a particular seismic phase, represented by the traveltime sensitivity kernel K ( t 0 , x s ) with respect to the earthquake location ( t 0 , x s ), is theoretically derived based on the adjoint method. Traveltime sensitivity kernel K ( t 0 , x s ) is formulated as a convolution between the forward and adjoint wavefields, which are calculated by numerically solving two acoustic wave equations. The advantage of this newly derived traveltime kernel is that it not only takes into account the earthquake–receiver geometry but also accurately honours the complexity of the velocity model. The earthquake location is obtained by solving a regularized least-squares problem. In 3-D realistic applications, it is computationally expensive to conduct full wave simulations. Therefore, we propose a 2.5-D approach which assumes the forward and adjoint wave simulations within a 2-D vertical plane passing through the earthquake and receiver. Various synthetic examples show the accuracy of this acoustic wave-equation-based earthquake location method. The accuracy and efficiency of the 2.5-D approach for 3-D earthquake location are further verified by its application to the 2004 Big Bear earthquake in Southern California.