ALBERT

Description: 〈span〉〈div〉Summary〈/div〉Seismic free oscillations, or normal modes, provide a convenient tool to calculate low-frequency seismograms in heterogeneous Earth models. A procedure called ‘full mode coupling’ allows the seismic response of the Earth to be computed. However, in order to be theoretically exact, such calculations must involve an infinite set of modes. In practice, only a finite subset of modes can be used, introducing an error into the seismograms. By systematically increasing the number of modes beyond the highest frequency of interest in the seismograms, we investigate the convergence of full-coupling calculations. As a rule-of-thumb, it is necessary to couple modes 1–2 mHz above the highest frequency of interest, although results depend upon the details of the Earth model. This is significantly higher than has previously been assumed. Observations of free oscillations also provide important constraints on the heterogeneous structure of the Earth. Historically, this inference problem has been addressed by the measurement and interpretation of splitting functions. These can be seen as secondary data extracted from low frequency seismograms. The measurement step necessitates the calculation of synthetic seismograms, but current implementations rely on approximations referred to as self- or group-coupling and do not use fully accurate seismograms. We therefore also investigate whether a systematic error might be present in currently published splitting functions. We find no evidence for any systematic bias, but published uncertainties must be doubled to properly account for the errors due to theoretical omissions and regularization in the measurement process. Correspondingly, uncertainties in results derived from splitting functions must also be increased. As is well known, density has only a weak signal in low-frequency seismograms. Our results suggest this signal is of similar scale to the true uncertainties associated with currently published splitting functions. Thus, it seems that great care must be taken in any attempt to robustly infer details of Earth's density structure using current splitting functions.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉In this paper, we propose a new wavelet-based 3-D inversion method for frequency-domain airborne electromagnetic (FDAEM) data. Instead of inverting the model in the space domain using a smoothing constraint, this new method recovers the model in the wavelet domain based on a sparsity constraint. In the wavelet domain, the model is represented by two types of coefficients, which contain both large- and fine-scale informations of the model, meaning the wavelet-domain inversion has inherent multiresolution. In order to accomplish a sparsity constraint, we minimize an L〈sub〉1〈/sub〉-norm measure in the wavelet domain that mostly gives a sparse solution. The final inversion system is solved by an iteratively reweighted least-squares method. We investigate different orders of Daubechies wavelets to accomplish our inversion algorithm, and test them on synthetic frequency-domain AEM data set. The results show that higher order wavelets having larger vanishing moments and regularity can deliver a more stable inversion process and give better local resolution, while the lower order wavelets are simpler and less smooth, and thus capable of recovering sharp discontinuities if the model is simple. At last, we test this new inversion algorithm on a frequency-domain helicopter EM (HEM) field data set acquired in Byneset, Norway. Wavelet-based 3-D inversion of HEM data is compared to L〈sub〉2〈/sub〉-norm-based 3-D inversion's result to further investigate the features of the new method.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉A monitoring method to grasp the spatio-temporal change in the interplate coupling in a subduction zone based on the spatial gradients of surface displacement rate fields is proposed. I estimated the spatio-temporal change in the interplate coupling along the plate boundary in northeastern (NE) Japan by applying the proposed method to the surface displacement rates based on global positioning system observations. The gradient of the surface velocities is calculated in each swath configured along the direction normal to the Japan Trench for time windows such as 0.5, 1, 2, 3 and 5 yr being shifted by one week during the period of 1997–2016. The gradient of the horizontal velocities is negative and has a large magnitude when the interplate coupling at the shallow part (less than approximately 50 km in depth) beneath the profile is strong, and the sign of the gradient of the vertical velocity is sensitive to the existence of the coupling at the deep part (greater than approximately 50 km in depth). The trench-parallel variation of the spatial gradients of a displacement rate field clearly corresponds to the trench-parallel variation of the amplitude of the interplate coupling on the plate interface, as well as the rupture areas of previous interplate earthquakes. Temporal changes in the trench-parallel variation of the spatial gradient of the displacement rate correspond to the strengthening or weakening of the interplate coupling. We can monitor the temporal change in the interplate coupling state by calculating the spatial gradients of the surface displacement rate field to some extent without performing inversion analyses with applying certain constraint conditions that sometimes cause over- and/or underestimation at areas of limited spatial resolution far from the observation network. The results of the calculation confirm known interplate events in the NE Japan subduction zone, such as the post-seismic slip of the 2003 M8.0 Tokachi-oki and 2005 M7.2 Miyagi-oki earthquakes and the recovery of the interplate coupling around the rupture area of the 1994 M7.6 Sanriku-Haruka-oki earthquake. The results also indicate the semi-periodic occurrence of slow slip events and the expansion of the area of slow slip events before the 2011 Tohoku-oki earthquake (M9.0) approaching the hypocentre of the Tohoku-oki earthquake.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We image the internal structure of the San Jacinto fault zone (SJFZ) in the trifurcation area southeast of Anza, California, with seismic records from dense linear and rectangular arrays. The examined data include recordings from more than 20 000 local earthquakes and nine teleseismic events. Automatic detection algorithms and visual inspection are used to identify 〈span〉P〈/span〉 and 〈span〉S〈/span〉 body waves, along with 〈span〉P〈/span〉- and 〈span〉S〈/span〉-types fault zone trapped waves (FZTW). The location at depth of the main branch of the SJFZ, the Clark fault, is identified from systematic waveform changes across lines of sensors within the dense rectangular array. Delay times of 〈span〉P〈/span〉 arrivals from teleseismic and local events indicate damage asymmetry across the fault, with higher damage to the NE, producing a local reversal of the velocity contrast in the shallow crust with respect to the large-scale structure. A portion of the damage zone between the main fault and a second mapped surface trace to the NE generates 〈span〉P〈/span〉- and 〈span〉S〈/span〉-types FZTW. Inversions of high-quality 〈span〉S〈/span〉-type FZTW indicate that the most likely parameters of the trapping structure are width of ∼70 m, 〈span〉S〈/span〉-wave velocity reduction of 60 per cent, 〈span〉Q〈/span〉 value of 60 and depth of ∼2 km. The local reversal of the shallow velocity contrast across the fault with respect to large-scale structure is consistent with preferred propagation of earthquake ruptures in the area to the NW.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The computational cost of quasi-〈span〉P〈/span〉 wave extrapolation depends on the complexity of the medium, and specifically the anisotropy. Our effective-model method splits the anisotropic dispersion relation into an isotropic background and a correction factor to handle this dependency. The correction term depends on the slope (measured using the gradient) of current wavefields and the anisotropy. As a result, the computational cost is independent of the nature of anisotropy, which makes the extrapolation efficient. A dynamic implementation of this approach decomposes the original pseudo-differential operator into a Laplacian, handled using the low-rank approximation of the spectral operator, plus an angular dependent correction factor applied in the space domain to correct for anisotropy. We analyse the role played by the correction factor and propose a new spherical decomposition of the dispersion relation. The proposed method provides accurate wavefields in phase and more balanced amplitudes than a previous spherical decomposition. Also, it is free of 〈span〉SV〈/span〉-wave artefacts. Applications to a simple homogeneous transverse isotropic medium with a vertical symmetry axis (VTI) and a modified Hess VTI model demonstrate the effectiveness of the approach. The Reverse Time Migration applied to a modified BP VTI model reveals that the anisotropic migration using the proposed modelling engine performs better than an isotropic migration.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Over the past 15 yr, numerical models of convection in Earth’s mantle have made a leap forward: they can now produce self-consistent plate-like behaviour at the surface together with deep mantle circulation. These digital tools provide a new window into the intimate connections between plate tectonics and mantle dynamics, and can therefore be used for tectonic predictions, in principle. This contribution explores this assumption. First, initial conditions at 30, 20, 10 and 0 Ma are generated by driving a convective flow with imposed plate velocities at the surface. We then compute instantaneous mantle flows in response to the guessed temperature fields without imposing any boundary conditions. Plate boundaries self-consistently emerge at correct locations with respect to reconstructions, except for small plates close to subduction zones. As already observed for other types of instantaneous flow calculations, the structure of the top boundary layer and upper-mantle slab is the dominant character that leads to accurate predictions of surface velocities. Perturbations of the rheological parameters have little impact on the resulting surface velocities. We then compute fully dynamic model evolution from 30 and 10 to 0 Ma, without imposing plate boundaries or plate velocities. Contrary to instantaneous calculations, errors in kinematic predictions are substantial, although the plate layout and kinematics in several areas remain consistent with the expectations for the Earth. For these calculations, varying the rheological parameters makes a difference for plate boundary evolution. Also, identified errors in initial conditions contribute to first-order kinematic errors. This experiment shows that the tectonic predictions of dynamic models over 10 My are highly sensitive to uncertainties of rheological parameters and initial temperature field in comparison to instantaneous flow calculations. Indeed, the initial conditions and the rheological parameters can be good enough for an accurate prediction of instantaneous flow, but not for a prediction after 10 My of evolution. Therefore, inverse methods (sequential or data assimilation methods) using short-term fully dynamic evolution that predict surface kinematics are promising tools for a better understanding of the state of the Earth’s mantle.〈/span〉

Description: Elastic reverse time migration (RTM) can yield accurate subsurface information (e.g. PP and PS reflectivity) by imaging the multicomponent seismic data. However, the existing RTM methods are still insufficient to provide satisfactory results because of the finite recording aperture, limited bandwidth and imperfect illumination. Besides, the P - and S -wave separation and the polarity reversal correction are indispensable in conventional elastic RTM. Here, we propose an iterative elastic least-squares RTM (LSRTM) method, in which the imaging accuracy is improved gradually with iteration. We first use the Born approximation to formulate the elastic de-migration operator, and employ the Lagrange multiplier method to derive the adjoint equations and gradients with respect to reflectivity. Then, an efficient inversion workflow (only four forward computations needed in each iteration) is introduced to update the reflectivity. Synthetic and field data examples reveal that the proposed LSRTM method can obtain higher-quality images than the conventional elastic RTM. We also analyse the influence of model parametrizations and misfit functions in elastic LSRTM. We observe that Lamé parameters, velocity and impedance parametrizations have similar and plausible migration results when the structures of different models are correlated. For an uncorrelated subsurface model, velocity and impedance parametrizations produce fewer artefacts caused by parameter crosstalk than the Lamé coefficient parametrization. Correlation- and convolution-type misfit functions are effective when amplitude errors are involved and the source wavelet is unknown, respectively. Finally, we discuss the dependence of elastic LSRTM on migration velocities and its antinoise ability. Imaging results determine that the new elastic LSRTM method performs well as long as the low-frequency components of migration velocities are correct. The quality of images of elastic LSRTM degrades with increasing noise.

Description: In general, the complex electrical resistivity in the subsurface is anisotropic. Despite this, algorithms for the tomographic inversion of complex resistivity data commonly assume isotropy, mainly due to the lack of anisotropic modelling and inversion schemes, potentially leading to artefacts in the inversion results in the presence of anisotropy. The development of an effective anisotropic complex resistivity inversion algorithm which utilizes the gradient information of some cost function benefits from understanding the characteristics of the problem's sensitivities, that is, the partial derivative of the impedance forward response with respect to the complex conductivities in the different spatial directions, as well as with respect to the different ratios of complex conductivities, that is, the different anisotropy ratios. We here derive expressions for these sensitivities and, based on a 2.5-D finite-element modelling algorithm, we compute and discuss sensitivity distributions as well as measurement response curves of typical surface and cross-borehole measurement configurations for 2-D subsurface anisotropic complex resistivity distributions. Depending on the electrode layout and measurement configuration, the sensitivity with respect to the conductivity in a particular direction shows a unique pattern, while for other directions sensitivity patterns are qualitatively similar. These sensitivity characteristics translate into important equivalences between impedance responses of local anisotropic and isotropic anomalies, for both magnitude and phase. Accordingly, with collinear surface arrays only the complex conductivity in the direction of the electrode layout can be unambiguously resolved, and with cross-borehole arrays only the conductivity in the vertical direction, provided an in-hole current injection is used. Nevertheless, anisotropy ratios involving these resolvable conductivity components are likewise detectable. The distinct shape of the measurement response curves, reflecting the distinct spatial patterns of the corresponding sensitivity distributions, suggest that optimized measurement configurations can be inferred for specific exploration questions involving electrical anisotropy and given electrode layouts. The gained insight into the characteristics of the sensitivity distributions of complex resistivity measurements in case of subsurface anisotropy should guide the implementation of effective anisotropic complex resistivity inversion schemes and lead to a routine use of such schemes in any resistivity and induced polarization surveys whenever subsurface electrical anisotropy could be encountered.

Description: The streaming potential phenomenon is an electrokinetic effect that occurs in porous media. It is characterized by an electrokinetic (EK) coefficient. The aim of this paper is to simulate the EK coefficient in unsaturated conditions using the Lattice Boltzmann method in a 2-D capillary channel. The multiphase flow is simulated with the model of Shan & Chen. The Poisson–Boltzmann equation is solved by implementing the model of Chai & Shi. The streaming potential response shows a non-monotonous behaviour due to the combination of the increase of charge density and decrease of flow velocity with decreasing water saturation. Using a potential of –20 mV at the air–water interface, an enhancement of a factor 5–30 of the EK coefficient, compared to the saturated state, can be observed due to the positive charge excess at this interface which is magnified by the fluid velocity away from the rock surface. This enhancement is correlated to the fractioning of the bubbles, and to the dynamic state of these bubbles, moving or entrapped in the crevices of the channel.

Description: Seismic noise measurements (ambient vibrations) have been increasingly used in rock slope stability assessment for both investigation and monitoring purposes. Recent studies made on gravitational hazard revealed significant spectral amplification at given frequencies and polarization of the wave-field in the direction of maximum rock slope displacement. Different properties (resonance frequencies, polarization and spectral ratio amplitudes) can be derived from the spectral analysis of the seismic noise to characterize unstable rock masses. The objective here is to identify the dynamic parameters that could be used to gain information on prone-to-fall rock columns’ geometry. To do so, the dynamic response of prone-to-fall columns to seismic noise has been studied on two different sites exhibiting cliff-like geometry. Dynamic parameters (main resonance frequency and spectral ratio amplitudes) that could characterize the column decoupling were extracted from seismic noise and their variations were studied taking into account the external environmental parameter fluctuations. Based on this analysis, a two-dimensional numerical model has been set up to assess the influence of the rear vertical fractures identified on both sites on the rock column motion response. Although a simple relation was found between spectral ratio amplitudes and the rock column slenderness, it turned out that the resonance frequency is more stable than the spectral ratio amplitudes to characterize this column decoupling, provided that the elastic properties of the column can be estimated. The study also revealed the effect of additional remote fractures on the dynamic parameters, which in turn could be used for detecting the presence of such discontinuities.

Description: In this study, we present a new synthesis of GPS velocities for tectonic deformation within the Tibetan Plateau and its surrounding areas, a combined data set of ~1854 GPS-derived horizontal velocity vectors. Assuming that crustal deformation is localized along major faults, a block modelling approach is employed to interpret the GPS velocity field. We construct a 30-element block model to describe present-day deformation in western China, with half of them located within the Tibetan Plateau, and the remainder located in its surrounding areas. We model the GPS velocities simultaneously for the effects of block rotations and elastic strain induced by the bounding faults. Our model yields a good fit to the GPS data with a mean residual of 1.08 mm a –1 compared to the mean uncertainty of 1.36 mm a –1 for each velocity component, indicating a good agreement between the predicted and observed velocities. The major strike-slip faults such as the Altyn Tagh, Xianshuihe, Kunlun and Haiyuan faults have relatively uniform slip rates in a range of 5–12 mm a –1 along most of their segments, and the estimated fault slip rates agree well with previous geologic and geodetic results. Blocks having significant residuals are located at the southern and southeastern Tibetan Plateau, suggesting complex tectonic settings and further refinement of accurate definition of block geometry in these regions.

Description: Surface-related multiples have been utilized in the reverse-time migration (RTM) procedures, and additional illumination for subsurface can be provided. Meanwhile, many cross-talks are generated from undesired interactions between forward- and backward-propagated seismic waves. In this paper, subsequent to analysing and categorizing these cross-talks, we propose RTM of first-order multiples to avoid most undesired interactions in RTM of all-order multiples, where only primaries are forward-propagated and crosscorrelated with the backward-propagated first-order multiples. With primaries and multiples separated during regular seismic data processing as the input data, first-order multiples can be obtained by a two-step scheme: (1) the dual-prediction of higher-order multiples; and (2) the adaptive subtraction of predicted higher-order multiples from all-order multiples within local offset-time windows. In numerical experiments, two synthetic and a marine field data sets are used, where different cross-talks generated by RTM of all-order multiples can be identified and the proposed RTM of first-order multiples can provide a very interpretable image with a few cross-talks.

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉We present the theory of parsimonious refraction interferometry and tomography where a densely populated refraction data set can be obtained from two reciprocal and several infill shot gathers. The assumptions are that the refraction arrivals are head waves, and a pair of reciprocal shot gathers and several infill shot gathers are recorded over the line of interest. Refraction traveltimes from these shot gathers are picked and spawned into 〈span style="font-style:italic;"〉O〈/span〉(〈span style="font-style:italic;"〉N〈/span〉〈sup〉2〈/sup〉) virtual refraction traveltimes generated by 〈span style="font-style:italic;"〉N〈/span〉 virtual sources, where 〈span style="font-style:italic;"〉N〈/span〉 is the number of geophones in the 2-D survey. The virtual traveltimes can be inverted to give the velocity tomogram. This enormous increase in the number of traveltime picks and associated rays, compared to the many fewer traveltimes from the reciprocal and infill shot gathers, allows for increased model resolution and a better condition number with the system of normal equations. A significant benefit is that the parsimonious survey and the associated traveltime picking can be an order-of-magnitude less time consuming than that for a standard refraction survey with a dense distribution of sources.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉We describe three achievements for a ground motion simulation. First, we propose a kinematic modelling in which rupture delay time is governed by an eikonal equation on a Riemannian manifold and develop a coupling method between the ground motion simulation and the eikonal solver. In general the rupture velocity distribution is not spatially uniform and the rupture propagation depends on a fault shape. So we derive the eikonal equation by considering the Riemannian metric of the fault surface and give a detailed discretization of its difference scheme to deal with a curved surface fault. Next, in order to explain the effect of spatially discontinuous non-uniformity of rupture velocity, we introduce an isochrones jumping intensity and obtain a new decomposed isochrones formula in general settings. It is known that the representation theorem with the Green's function can be rewritten into an expression with a contour integral by the isochrones theory. The new formula says that the known isochrones formula for ground velocity can be decomposed into a trend component and a disturbance component. The disturbance component consists of the isochrones jumping intensity. Finally, by applying our ground motion simulation coupled with the eikonal solver and the decomposed isochrones formula, we investigate some relations between the non-uniformity of the rupture velocity and pulse-like disturbance of the ground motion velocity. Our simulations show that the disturbance of velocity waveform corresponds with that of rate of change of isochrones band area. It turns out that the pulse-like disturbance of velocity waveform occurs when isochrones move across the region where rupture velocity varies discontinuously. Thus we can explain that the pulse-like disturbance of the ground motion velocity occurs when the isochrones jumping intensity has nonzero value. Moreover, as another example of application of our simulation and formula, we show a distinctive dependence of peak ground velocity upon parameters such as the rupture velocity and the distance between a fault and an observer.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉Over the last decades, electromagnetic methods have become an accepted tool for a wide range of geophysical exploration purposes and nowadays even for monitoring. Application to hydrocarbon monitoring, for example for enhanced oil recovery, is hampered by steel-cased wells, which typically exist in large numbers in producing oil fields and which distort electromagnetic fields in the subsurface. Steel casings have complex geometries as they are very thin but vertically extended; moreover, the conductivity contrast of steel to natural materials is in the range of six orders of magnitude. It is therefore computationally prohibitively costly to include such structures directly into the modelling grid, even for finite element methods. To tackle the problem we developed a method to describe steel-cased wells as series of substitute dipole sources, which effectively interact with the primary field. The new approach cannot only handle a single steel-cased well, but also an arbitrary number, and their interaction with each other. We illustrate the metal casing effect with synthetic 3-D modelling of land-based controlled source electromagnetic data. Steel casings distort electromagnetic fields even for large borehole-transmitter distances above 2 km. The effect depends not only on the distance between casing and transmitter, but also on the orientation of the transmitter to the borehole. Finally, we demonstrate how the presence of steel-cased wells can be exploited to increase the sensitivity and enhance resolution in the target region. Our results show that it is at least advisable to consider the distribution of steel-cased wells already at the planning phase of a controlled source electromagnetic field campaign.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉The volatile content in magmas is fundamental for the triggering and style of volcanic eruptions. Carbon dioxide, the second most abundant volatile component in magmas after H〈sub〉2〈/sub〉O, is the first to reach saturation upon ascent and depressurization. We investigate experimentally CO〈sub〉2〈/sub〉-bubble nucleation in trachybasalt and trachyte melts at high temperature and high pressure (〈span style="font-style:italic;"〉HT〈/span〉 and 〈span style="font-style:italic;"〉HP〈/span〉) through wetting-angle measurements on different (sialic, mafic or oxide) phenocryst phases. The presence of crystals lowers the supersaturation required for CO〈sub〉2〈/sub〉-bubble nucleation up to 37 per cent (heterogeneous nucleation, 〈span style="font-style:italic;"〉HeN〈/span〉), with a minor role of mineral chemistry. Different from H〈sub〉2〈/sub〉O-rich systems, feldspar crystals are effective in reducing required supersaturation for bubble nucleation. Our data suggest that leucite, the dominant 〈span style="font-style:italic;"〉liquidus〈/span〉 phase in ultrapotassic systems at shallow depth (i.e. 〈100 MPa), facilitates late-stage, extensive magma vesiculation through CO〈sub〉2〈/sub〉〈span style="font-style:italic;"〉HeN〈/span〉, which may explain the shifting of CO〈sub〉2〈/sub〉-rich eruptive systems towards an apparently anomalous explosive behaviour.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The contribution of the GOCE gravity gradients to regional gravity field solutions is investigated in this study. We employ radial basis functions to recover the gravity field on regional scales over Amazon and Himalayas as our test regions. In the first step, four individual solutions based on the more accurate gravity gradient components 〈span style="font-style:italic;"〉Txx〈/span〉, 〈span style="font-style:italic;"〉Tyy〈/span〉, 〈span style="font-style:italic;"〉Tzz〈/span〉 and 〈span style="font-style:italic;"〉Txz〈/span〉 are derived. The 〈span style="font-style:italic;"〉Tzz〈/span〉 component gives better solution than the other single-component solutions despite the less accuracy of 〈span style="font-style:italic;"〉Tzz〈/span〉 compared to 〈span style="font-style:italic;"〉Txx〈/span〉 and 〈span style="font-style:italic;"〉Tyy〈/span〉. Furthermore, we determine five more solutions based on several selected combinations of the gravity gradient components including a combined solution using the four gradient components. The 〈span style="font-style:italic;"〉Tzz〈/span〉 and 〈span style="font-style:italic;"〉Tyy〈/span〉 components are shown to be the main contributors in all combined solutions whereas the 〈span style="font-style:italic;"〉Txz〈/span〉 adds the least value to the regional gravity solutions. We also investigate the contribution of the regularization term. We show that the contribution of the regularization significantly decreases as more gravity gradients are included. For the solution using all gravity gradients, regularization term contributes to about 5 per cent of the total solution. Finally, we demonstrate that in our test areas, regional gravity modelling based on GOCE data provide more reliable gravity signal in medium wavelengths as compared to pre-GOCE global gravity field models such as the EGM2008.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Defining variations in the behaviour of the geomagnetic field through geological time is critical to understanding the dynamics of Earth's core and its response to mantle convection and planetary evolution. Furthermore, the question of whether the axial dipole dominance of the recent palaeomagnetic field persists through the whole of Earth's history is fundamental to determining the reliability of palaeogeographic reconstructions and the efficacy of the magnetosphere in shielding Earth from solar wind radiation. Previous palaeomagnetic directional studies have suggested that the palaeofield had a complex configuration in the Devonian period (419–359 Ma). Here we present new high-quality palaeointensity determinations from rocks aged between 408 and 375 Ma from the Minusa Basin (southern Siberia), and the Kola Peninsula that enable the first reliable investigation of the strength of the field during this enigmatic period. Palaeointensity experiments were performed using the thermal Thellier, microwave Thellier and Wilson methods on 165 specimens from 25 sites. Six out of eight successful sites from the Minusa Basin and all four successful sites from the Kola Peninsula produced extremely low palaeointensities (〈10 μT). These findings challenge the uniformitarian view of the palaeomagnetic field: field intensities of nearly an order of magnitude lower than Neogene values (except during relatively rare geomagnetic excursions and reversals) together with the widespread appearance of strange directions found in the Devonian suggest that the Earth's field during this time may have had a dominantly multipolar geometry. A persistent, low intensity multipolar magnetic field and associated diminished magnetosphere would increase the impact of solar particles on the Earth's magnetosphere, ionosphere and atmosphere with potential major implications for Earth's climate and biosphere.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Glacial Isostatic Adjustment (GIA) models commonly assume a mantle with a viscoelastic Maxwell rheology and a fixed ice history model. Here, we use a Bayesian Monte Carlo approach with a Markov chain formalism to invert the global GIA signal simultaneously for the mechanical properties of the mantle and the volumes of the ice sheets, using as starting ice models two previously published ice histories. Two stress relaxing rheologies are considered: Burgers and Maxwell linear viscoelasticities. A total of 5720 global palaeo sea level records are used, covering the last 35 kyr. Our goal is not only to seek the model best fitting this data set, but also to determine and display the range of possible solutions with their respective probability of explaining the data. In all cases, our 〈span style="font-style:italic;"〉a posteriori〈/span〉 probability maps exhibit the classic character of solutions for GIA-determined mantle viscosity with two distinct peaks. What is new in our treatment is the presence of the bi-viscous Burgers rheology and the fact that we invert rheology jointly with ice history, in combination with the greatly expanded palaeo sea level records. The solutions tend to be characterized by an upper-mantle viscosity of around 5 × 10〈sup〉20〈/sup〉 Pa s with one preferred lower-mantle viscosities at 3 × 10〈sup〉21〈/sup〉 Pa s and the other more than 2 × 10〈sup〉22〈/sup〉 Pa s, a rather classical pairing. Best-fitting models depend upon the starting ice history and the stress relaxing law. A first peak (P1) has the highest probability only in the case with a Maxwell rheology and ice history based on ICE-5G, while the second peak (P2) is favoured for ANU-based ice history or Burgers stress relaxation. The latter solution also may satisfy lower-mantle viscosity inferences from long-term geodynamics and gravity gradient anomalies over Laurentia. P2 is also consistent with large Laurentian and Fennoscandian ice-sheet volumes at the Last Glacial Maximum (LGM) and smaller LGM Antarctic ice volume than in either ICE-5G or ANU. Exploration of a bi-viscous linear relaxing rheology in GIA now seems logical due to a new set of requirements to satisfy observations of transient post-seismic flow seen so ubiquitously in space gravimetry and other global geodetic data.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Earthquakes mainly occur in crust or mantle that is below a critical temperature for the tectonic strain-rate, $\dot{e}_t$, such that stress builds up to the breaking point before it can relax due to creep. Then long-range stress correlation gives rise to power law seismicity including large events. The limiting temperature depends on pressure, which is taken into account by finding a critical homologous temperature 〈span style="font-style:italic;"〉T〈/span〉〈sub〉Hc〈/sub〉 = 〈span style="font-style:italic;"〉T〈/span〉/〈span style="font-style:italic;"〉T〈/span〉〈sub〉M〈/sub〉 above which earthquakes are rarely observed (where 〈span style="font-style:italic;"〉T〈/span〉, 〈span style="font-style:italic;"〉T〈/span〉〈sub〉M〈/sub〉 are temperature and average melting temperature of constituent minerals). We find that 〈span style="font-style:italic;"〉T〈/span〉〈sub〉Hc〈/sub〉 for ocean plates is ∼0.55. For California earthquakes, it is also close to 0.55. The uppermost mantle layer of oceanic plates of thickness ∼50 km is composed of harzburgite and depleted peridotite from which basalt has been removed to form ocean crust. Thus it has a higher melting temperature than the peridotite of the surrounding mantle, or the lower halves of plates. Thicknesses of seismicity in deep subduction zones, determined from 2-D polynomial fits to a relocated catalogue, are ∼50 km, which suggests that the earthquake channel is confined to this layer. We construct models to find homologous temperatures in slabs, and find that seismicity thicknesses are also, on average, confined to 〈span style="font-style:italic;"〉T〈/span〉〈sub〉H〈/sub〉 ≤ 0.55 ± 0.05. The associated rheology is compared with that obtained from flexure models of ocean lithosphere. The brittle-ductile transition occurs where viscosity drops from high values in the cold cores of slabs to values of 10〈sup〉22〈/sup〉–10〈sup〉23〈/sup〉 Pa s, that is, where creep strain-rates become comparable to tectonic rates. The cut-off for deep earthquakes is not sharp. However they appear unlikely to occur if homologous temperature is high 〈span style="font-style:italic;"〉T〈/span〉〈sub〉H〈/sub〉 〉 0.55. Exceptions to the rule are anomalously deep earthquakes such as those beneath the Iceland and the Hawaiian hotspots, and the Newport Inglewood Fault. These are smaller events with short-range stress correlation, and can be explained if strain-rates are two to three orders of magnitude higher than those associated with earthquakes located where 〈span style="font-style:italic;"〉T〈/span〉〈sub〉H〈/sub〉 ≤ 0.55. We conclude that the brittle-ductile transition corresponds to the transition from long-range (regional) to short-range (localized on asperities) stress correlation.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We introduce a wavefield gradiometry technique to estimate both isotropic and anisotropic local medium characteristics from short recordings of seismic signals by inverting a wave equation. The method exploits the information in the spatial gradients of a seismic wavefield that are calculated using dense deployments of seismic arrays. The application of the method uses the surface wave energy in the ambient seismic field. To estimate isotropic and anisotropic medium properties we invert an elliptically anisotropic wave equation. The spatial derivatives of the recorded wavefield are evaluated by calculating finite differences over nearby recordings, which introduces a systematic anisotropic error. A two-step approach corrects this error: finite difference stencils are first calibrated, then the output of the wave-equation inversion is corrected using the linearized impulse response to the inverted velocity anomaly. We test the procedure on ambient seismic noise recorded in a large and dense ocean bottom cable array installed over Ekofisk field. The estimated azimuthal anisotropy forms a circular geometry around the production-induced subsidence bowl. This conforms with results from studies employing controlled sources, and with interferometry correlating long records of seismic noise. Yet in this example, the results were obtained using only a few minutes of ambient seismic noise.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Reverberations of teleseismic compressional (〈span style="font-style:italic;"〉P〈/span〉-) waves within a glacier or ice sheet may mask signals associated with crustal structure beneath the ice. We remove the signal associated with the ice from teleseismic 〈span style="font-style:italic;"〉P〈/span〉-waves using a wavefield downward continuation and decomposition technique that depends on known ice layer properties such as ice thickness, velocity, and attenuation. We test the method using data from nine stations in Antarctica and one station in Greenland. We deconvolve the downward-continued seismic wave vectors to create 〈span style="font-style:italic;"〉P〈/span〉-wave receiver functions that minimize the ice-layer reverberations in order to better measure signals from deeper structures. The subsurface 〈span style="font-style:italic;"〉P〈/span〉-wave receiver functions have similar sensitivities to crustal structure as those calculated from stations installed on bedrock. Synthetic experiments indicate subsurface 〈span style="font-style:italic;"〉P〈/span〉-wave receiver functions can constrain crustal structure more tightly than surface 〈span style="font-style:italic;"〉P〈/span〉-wave receiver functions when ice layer properties are known. We model the subsurface 〈span style="font-style:italic;"〉P〈/span〉-wave receiver functions using a Markov chain Monte Carlo inversion and constrain the product of crustal thickness and the column-average crustal-slowness beneath the stations. Our subglacial shear speed and thickness estimates are consistent with previous investigations at most stations. At station SUMG in south-central Greenland, our results suggest a thicker crust than from previous estimates.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Usually, when inverting geodetic data to estimate the slip distributions on a fault, the area is made large enough to more than cover the rupture zone, with regularization producing regions of large slip with very small slip over the rest of the surface. We have developed a new inverse method which assumes that nonzero slip is confined to a rectangular region, and which jointly estimates, using Bayesian methods, the boundaries of this region as well as the slip distribution within it, using a smoothing parameter also determined as part of the inversion. Synthetic tests show that our method can successfully image deeper slip regions not resolved by previous methods, and does not produce spurious regions of nonzero slip. We apply our method to coseismic displacements measured by GPS for the 2009 L’Aquila earthquake, first determining the orientation of the fault assuming a simplified model with uniform slip, and then determining probability density functions for the location, length, and width of the rupture area and for the slip distribution. The standard deviation of slip is about 10 cm and describes a normal-faulting earthquake with a maximum slip of 88 ± 11 cm and seismic moment of $3.32_{-0.29}^{+0.30}\times 10^{18}$ N m.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉High-resolution models of seismic velocity variations constructed using body-wave tomography inform the study of the origin, fate and thermochemical state of mantle domains. In order to reliably relate these variations to material properties including temperature, composition and volatile content, we must accurately retrieve both the patterns and amplitudes of variations and quantify the uncertainty associated with the estimates of each. For these reasons, we image the mantle beneath North America with 〈span style="font-style:italic;"〉P〈/span〉-wave traveltimes from USArray using a novel method for 3-D probabilistic body-wave tomography. The method uses a Transdimensional Hierarchical Bayesian framework with a reversible-jump Markov Chain Monte Carlo algorithm in order to generate an ensemble of possible velocity models. We analyse this ensemble solution to obtain the posterior probability distribution of velocities, thereby yielding error bars and enabling rigorous hypothesis testing. Overall, we determine that the average uncertainty (1σ) of compressional wave velocity estimates beneath North America is ∼0.25 per cent 〈span style="font-style:italic;"〉dVP〈/span〉/〈span style="font-style:italic;"〉VP〈/span〉, increasing with proximity to complex structure and decreasing with depth. The addition of USArray data reduces the uncertainty beneath the Eastern US by over 50 per cent in the upper mantle and 25–40 per cent below the transition zone and ∼30 per cent throughout the mantle beneath the Western US. In the absence of damping and smoothing, we recover amplitudes of variations 10–80 per cent higher than a standard inversion approach. Accounting for differences in data coverage, we infer that the length scale of heterogeneity is ∼50 per cent longer at shallow depths beneath the continental platform than beneath tectonically active regions. We illustrate the model trade-off analysis for the Cascadia slab and the New Madrid Seismic Zone, where we find that smearing due to the limitations of the illumination is relatively minor.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We demonstrate with synthetic and field data that with sufficiently dense sampling wave-equation-based methods such as reverse time migration (RTM), implicitly forming array receiver functions (ARFs), perform better resolution wise than migration of common conversion point (CCP) stacks of traditional receiver functions. However, even with modern array deployments the sampling requirement is typically not met for teleseismic (earthquake) data. To enable RTM imaging with sparsely (and irregularly) sampled wavefields at the surface, we use an intermediate reconstruction based on sparsity promoting optimization using a curvelet (or wave packet) representation of the data, as an important and necessary pre-processing step. To suppress artefacts, the curvelet coefficients are constrained to represent the range of known directions present in the data. We show that our proposed pre-processing procedure (which may be viewed as generating ‘missing’ traces) can produce artefact-free data for RTM even if only 20 per cent of necessary data are available in the original data set. With synthetic data, we also demonstrate that if the sampling criteria is not met, CCP can produce results that are superior over wave-equation methods such as RTM. As a proof-of-concept with field data, we image the structure of the crust beneath the Himalayas with passive-source RTM of teleseismic data from Hi-CLIMB project. For Hi-CLIMB data, the CCP and RTM results are similar because sampling is still too sparse for RTM and the structure is simple enough for successful CCP. Both results are improved by wavefield regularization and reveal that the Moho is continuous beneath most of the array, and not fragmented as suggested by some earlier studies.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Elastic full waveform inversion (EFWI) aims to reduce the misfit between recorded and modelled multicomponent seismic data for deducing a detailed model of elastic parameters in the subsurface. Because the explicit computation and inversion of the Hessian matrix is extremely resource intensive, a gradient-based (rather than Hessian-based) minimization is generally applied for large-scale applications. However, the multiparameter trade-off effects cause cross-talks in the computed gradients and thus severely affect the convergence and the quality of the inverted model. Recently, preconditioning the gradients based on elastic wave mode decomposition has been suggested for mitigating the parameter trade-offs in the EFWI process. In this paper, we propose a mode decomposition (MD)-based EFWI approach in which the preconditioned gradients are obtained through the cross-correlation of the forward and decomposed adjoint wavefields in the time domain. Based on the decomposed Frechét derivatives, we explain the mechanism of this approach through analyses of Hessian and resolution matrices and comparison with the Gauss–Newton gradients. Numerical examples of a simple fluid-saturated model and the Marmousi-II model demonstrate that the MD-based preconditioned conjugate-gradient approach can mitigate the trade-off between the 〈span style="font-style:italic;"〉P〈/span〉- and 〈span style="font-style:italic;"〉S〈/span〉-wave velocities and achieve fast convergence without any Hessian-involved calculations.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We use frequency domain methods usually applied to volcanic tremor to analyse ground based seismic recordings of a helicopter. We preclude misinterpretations of tremor sources and show alternative applications of our seismological methods. On a volcano, the seismic source can consist of repeating, closely spaced, small earthquakes. Interestingly, similar signals are generated by helicopters due to repeating pressure pulses from the rotor blades. In both cases the seismic signals are continuous and referred to as tremor. As frequency gliding is in this case merely caused by the Doppler effect, not a change in the source, we can use its shape to deduce properties of the helicopter and its flight path. We show in this analysis that the number of rotor blades, rotor revolutions per minute, helicopter speed, flight direction, altitude and location can be deduced from seismometer recordings. Access to GPS determined flight path data from the helicopter offers us a robust way to test our location method.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉An analytical solution is derived concerning the linear run-up for any given initial wave generated over a sloping bathymetry. Due to the simplicity of the linear formulation, complex transformations are unnecessary, hence the shoreline motion is directly obtained in terms of the initial wave. This result supports not only maximum run-up invariance between linear and nonlinear theories but also the time evolution of shoreline motion and velocity, exhibiting good agreement with the nonlinear theory. The present formulation also allows quantifying the shoreline motion numerically from a customized initial waveform, including non-smooth functions. This is useful for numerical tests, laboratory experiments or realistic cases in which the initial disturbance might be retrieved from seismic data rather than using a theoretical model. It is also shown that the run-up calculation for the real case studied is consistent with the field observations.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉In studies of the magnetic properties of soils, the frequency-dependent magnetic susceptibility percentage (〈span style="font-style:italic;"〉χ〈/span〉〈sub〉FD〈/sub〉%) is often used for the identification of ultrafine magnetically superparamagnetic/stable single-domain (SP/SSD) particles. This parameter is commonly used as an indicator for increased pedogenesis. In strongly magnetic soils, the SP/SSD magnetic signal (mostly bio-pedogenic) may be masked by lithological signals; making pedogenesis hard to detect. In this study, we compare results for the detection of ultrafine SP/SSD magnetic particles in andic soils using two instruments: a Bartington MS2B dual-frequency meter and an AGICO Kappabridge MFK1-FA. In particular, the study focuses on the effect of pedogenesis by investigating the relationship between specific soil magnetic and chemical properties (soil organic carbon and pH〈sub〉H2O〈/sub〉). The values of 〈span style="font-style:italic;"〉χ〈/span〉〈sub〉FD〈/sub〉% obtained with the MS2B varied from 2.4 to 5.9 per cent, and mass-specific magnetic susceptibility (〈span style="font-style:italic;"〉χ〈/span〉〈sub〉LF〈/sub〉) from 283 to 1688 × 10〈sup〉−8〈/sup〉 m〈sup〉3〈/sup〉 kg〈sup〉−1〈/sup〉, while values of 〈span style="font-style:italic;"〉χ〈/span〉〈sub〉FD〈/sub〉% and 〈span style="font-style:italic;"〉χ〈/span〉〈sub〉LF〈/sub〉 obtained with the MFK1-FA varied from 2.7 to 8.2 per cent and from 299 to 1859 × 10〈sup〉−8〈/sup〉 m〈sup〉3〈/sup〉 kg〈sup〉−1〈/sup〉, respectively. Our results suggest that the detection of the SP/SSD magnetic fraction can be accomplished by comparing relative trends of 〈span style="font-style:italic;"〉χ〈/span〉〈sub〉FD〈/sub〉% along the soil profile. Moreover, the discrimination between bio-pedogenic and lithogenic magnetic contributions in the SP/SSD fraction is possible by comparing the 〈span style="font-style:italic;"〉χ〈/span〉〈sub〉FD〈/sub〉% and 〈span style="font-style:italic;"〉χ〈/span〉〈sub〉LF〈/sub〉 data determined in the fine earth (〈2 mm) and the coarse fraction (4–10 mm) samples down the soil profile.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We present and use a phase coherence approach to identify seismic signals that have similar path effects but different source time functions: co-located earthquakes and tremor. The method used is a phase coherence-based implementation of empirical matched field processing, modified to suit tremor analysis. It works by comparing the frequency-domain phases of waveforms generated by two sources recorded at multiple stations. We first cross-correlate the records of the two sources at a single station. If the sources are co-located, this cross-correlation eliminates the phases of the Green’s function. It leaves the relative phases of the source time functions, which should be the same across all stations so long as the spatial extent of the sources are small compared with the seismic wavelength. We therefore search for cross-correlation phases that are consistent across stations as an indication of co-located sources. We also introduce a method to obtain relative locations between the two sources, based on back-projection of interstation phase coherence. We apply this technique to analyse two tremor-like signals that are thought to be composed of a number of earthquakes. First, we analyse a 20 s long seismic precursor to a 〈span style="font-style:italic;"〉M〈/span〉 3.9 earthquake in central Alaska. The analysis locates the precursor to within 2 km of the mainshock, and it identifies several bursts of energy—potentially foreshocks or groups of foreshocks—within the precursor. Second, we examine several minutes of volcanic tremor prior to an eruption at Redoubt Volcano. We confirm that the tremor source is located close to repeating earthquakes identified earlier in the tremor sequence. The amplitude of the tremor diminishes about 30 s before the eruption, but the phase coherence results suggest that the tremor may persist at some level through this final interval.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We present a new method to improve the convergence of the well-known Parker's formula for the modelling of gravity and magnetic fields caused by sources with complex topography. In the original Parker's formula, two approximations are made, which may cause considerable numerical errors and instabilities: (1) the approximation of the forward and inverse continuous Fourier transforms using their discrete counterparts, the forward and inverse Fast Fourier Transform (FFT) algorithms; (2) the approximation of the exponential function with its Taylor series expansion. In a previous paper of ours, we have made an effort addressing the first problem by applying the Gauss-FFT method instead of the standard FFT algorithm. The new Gauss-FFT based method shows improved numerical efficiency and agrees well with space-domain analytical or hybrid analytical-numerical algorithms. However, even under the simplifying assumption of a calculation surface being a level plane above all topographic sources, the method may still fail or become inaccurate under certain circumstances. When the peaks of the topography approach the observation surface too closely, the number of terms of the Taylor series expansion needed to reach a suitable precision becomes large and slows the calculation. We show in this paper that this problem is caused by the second approximation mentioned above, and it is due to the convergence property of the Taylor series expansion that the algorithm becomes inaccurate for certain topographic models with large amplitudes. Based on this observation, we present a modified Parker's method using low rank approximation of the exponential function in virtue of the Chebfun software system. In this way, the optimal rate of convergence is achieved. Some pre-computation is needed but will not cause significant computational overheads. Synthetic and real model tests show that the method now works well for almost any practical topographic model, provided that the assumption, that the entire topographic mass lies below the observation surface, is met.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Seismic velocities in the sediments containing gas hydrates show a marked increase compared to the background, intuitively implying that there would be a reduction in the seismic attenuation (Q〈sup〉−1〈/sup〉) of such sediments. However, the attenuation measurements carried out in the sonic frequency range from various gas hydrate provinces show that there is a notable increase in seismic attenuation within the gas hydrate layers and the results obtained are counter-intuitive. In this work we try to compare the attenuation derived by applying frequency shift method to the multichannel seismic (MCS) and sonic data sets at the same location in the Krishna-Godavari (KG) basin. The role of complex geology in attenuating the seismic signal is also studied by generating synthetic seismic data for different geological models and by computing the corresponding attenuation. It has been found that the Q〈sup〉−1〈/sup〉 obtained from the field seismic data compares well with the Q〈sup〉−1〈/sup〉 computed for a simple layered geological model. The results indicate that the Q〈sup〉−1〈/sup〉 (0.0029) from field seismic data is ∼4.3 times lower than Q〈sup〉−1〈/sup〉 (0.0123–0.0125) obtained from the sonic data, implying that the thickness of the gas hydrate layer in the KG basin is not sufficient enough to average the bulk properties. Our results also indicate that there could be substantial contribution of the pore scale interaction to the observed attenuation.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Palaeomagnetic and rock magnetic studies of impact-related rocks can provide important constraints for deciphering geophysical records from suspected impact structures, their geochronology, and, in the case of very large impacts, their effect on the Earth as a whole. However, the palaeomagnetic record in impact-related rocks may be ambiguous because of the uncertain origin of their natural remanent magnetization (NRM). Towards this end, we carried out a comprehensive rock magnetic and mineralogical study of tagamites (impact melts) from the Jänisjärvi astrobleme, Russian Karelia. Chemical composition of magnetic minerals and non-magnetic matrix was evaluated by scanning electron microscopy (SEM) and X-ray analysis. Magnetic minerals were identified using thermomagnetic analysis at high and low temperatures, whereas their domain state was evaluated from hysteresis measurements and magnetic force microscopy. Jänisjärvi tagamites appear to belong to two essentially different types arising from the differences in the impact melt crystallization conditions. Type I tagamites were likely formed by an extremely rapid cooling of a superhot melt with initial temperatures well above 2000 °C. Type II tagamites originate from cooler and more iron-enriched melt. Common to the two types is that they both contain a substantial amount of fine inclusions in silicate matrix tens of nanometres to few micrometres in size, which appear to be a major, in some cases dominant, magnetic mineral carrying a significant part of rocks NRM. Structurally, these inclusions are heterogeneous objects consisting of two phases showing both chemical and magnetic contrast.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Previous work has demonstrated that geoelectrical measurements, acquired either along the Earth’s surface or in boreholes, can be sensitive to the presence of fractures. However, a lack of numerical approaches that are well suited to modelling electric current flow in fractured media prevents us from systematically exploring the links between geoelectrical measurements and fractured rock properties. To address this issue, we present a highly computationally efficient methodology for the numerical simulation of geoelectrical data in 2.5-D in complex fractured domains. Our approach is based upon a discrete-dual-porosity formulation, whereby the fractures and rock matrix are treated separately and coupled through the exchange of electric current between them. We first validate our methodology against standard analytical and finite-element solutions. Subsequent use of the approach to simulate geoelectrical data for a variety of different fracture configurations demonstrates the sensitivity of these data to important parameters such as the fracture density, depth, and orientation.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We study the seismic source of the 2015 (〈span style="font-style:italic;"〉M〈/span〉〈sub〉w〈/sub〉 6.7) Jujuy, Argentina intermediate depth earthquake. We first constrain the fault plane by using a teleseismic inversion and by determining the aftershock distribution. Then, we perform kinematic and dynamic inversions to retrieve the parameters that control the rupture process, using data at regional distances, and modelling the source as an elliptical patch. Best models suggest a subshear rupture propagation with a duration of ∼5 s. Results from the dynamic modelling suggest a stress drop of 11.87 MPa and a fracture energy rate of 2.95 MJ m〈sup〉−2〈/sup〉, which are slightly less but of the same order as those of other events of similar size. Finally, we perform a Monte-Carlo inversion to explore the behaviour of the frictional parameters in the solution space, and then we compare our results with other intraslab events. We find that the 〈span style="font-style:italic;"〉κ〈/span〉 parameter (ratio between strain energy and fracture energy) and the relation between seismic moment and stress drop are similar for all the considered events.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉The interferometric synthetic aperture radar (InSAR) data from the Japan Aerospace Exploration Agency ALOS-2 satellite show possible deformation associated with the 2016 January 6 North Korean nuclear test whereas the European Space Agency Sentinel-1A data are decorrelated. This is the first time that deformation related to a nuclear test has been measured since 1992. Here, I present two interpretations of the observed deformation: First, the deformation can be explained by a triggered landslide on the western slope of Mt Mantap, with a displacement of up to 10 cm across a patch of 1 km〈sup〉2〈/sup〉. Second, the observation may be from uplift created by the nuclear explosion. In the second interpretation, the location, depth and cavity size can be estimated from a topography-corrected homogenous half-space model (Mogi). The preferred location of the 2016 January 6 event is 41.2993°N 129.0715°E, with an uncertainty of 100 m. The estimated depth is 420–700 m and the cavity radius is 23–27 m. Based on empirical data and the assumption of granite as the host rock, the yield is estimated to be 11.6–24.4 kilotons of TNT, which is consistent with previous results based on seismic data. With these two interpretations, I demonstrate that InSAR data provide an independent tool to locate and estimate source characteristics of nuclear tests in North Korea. The ambiguity of interpretation is mainly due to the limited InSAR data acquisition. Future frequent data collection by current and upcoming InSAR satellites will allow full use of InSAR for nuclear monitoring and characterization in North Korea and around the world.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We present azimuthally anisotropic Rayleigh group velocity models from 8 to 35 s both offshore and onshore of the South Island of New Zealand. We use MOANA (Marine Observations of Anisotropy Near Aotearoa) broad-band ocean seismic data in combination with on land data from the New Zealand National Seismography Network to investigate the seismic structure of the flanks of the Australian–Pacific plate boundary. At 8 s, we observe low offshore group velocities best explained by the influence of the water layer and thick water-laden sediments. At long periods (20–30 s), group velocities are lower on the South Island relative to its offshore flanks, due to thickened crust beneath the island, with the lowest velocities primarily beneath the Southern Alps. Group velocity azimuthal anisotropy fast directions near the Alpine Fault align with the direction of relative plate motion between the Australian and Pacific plates. In the southern portion of the island, fast directions rotate anticlockwise, likely in response to a decrease in dextral shearing away from the plate boundary. Azimuthal anisotropy fast directions align with absolute plate motion offshore on the Pacific plate. Based on the depth sensitivity of our observations, we suggest diffuse deformation occurs throughout the crust. Our observations match trends in previous Pn anisotropy and SKS shear wave splitting observations, and therefore suggest a consistent pattern of distributed deformation throughout the lithosphere.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Significant lateral and depth variations of the inner core's properties, such as the large-scale hemispherical pattern, have been confirmed by a variety of seismological observations. However it is still unclear which dynamic processes in the core are responsible for these variations. Small-scale volumetric heterogeneity has been detected in the top layer of the inner core by PKiKP coda observations. Studies of these small-scale heterogeneities can provide critical information, such as the degree of alignment of iron crystals, the presence of possible partial melt and the grain size of iron crystals, all of which can be used to constrain the dynamic processes of the inner core. However, most previous observations sampled the inner core beneath the Pacific Ocean and Asia, often in the inner core's ‘eastern hemisphere’. We use seismic stations in the North America, including the Earthscope Transportable Array, to look at PKiKP and its coda waves. We find 21 events with clear signals. In agreement with previous studies, inner core scattering (ICS), resulting in clear PKiKP coda, is found at epicentral distances of 60°–95°. However, the ICS we observe in these 21 western hemisphere events is weaker than previously reported for the eastern hemisphere. Comparing our observations with numerical simulations, we conclude that this relatively weak ICS indicates small-scale heterogeneity in at least the top layer of the inner core beneath Central America. Combining our clear observations with previous studies suggests either a hemispherical difference, or a regional variation, of small-scale heterogeneity in the inner core.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The Qinling Orogen Belt (QOB), Weihe Graben (WG) and the southern margin of the Ordos Block (SOB), lying on the central portion of China, had been involved into the amalgamation of China (Asian continent) through subduction, collision to exhumation processes. The Moho fabrics beneath this region, recorded part of the evolution. Therefore, its thickness and internal structure may provide significant knowledge and contribute to the understanding the intracontinental deformation of central China. In this paper, in order to place constrain on the nature beneath the study area, nine large dynamite shots (the charge ≥500 kg) used to infer the internal structure and characteristics of the crustal boundary. We analyse the specific characteristics of the Moho reflection, the amplitude decay curves in near vertical zone and generate a single-fold profile; in addition, it also address the internal structure and discuss its implication. The Moho is approximately at the depth of 39 km beneath the North Qinling Orogen (NQB) and the WG, and at the depth of 42 km beneath the SOB. The Moho shows a subtle uplift and the crust is thin under the NQB. The north-dipping reflectors between the lower crust and the uppermost mantle extend to the middle of the WG, and the south-dipping reflectors in the lower crust of the NQB are truncated by the Moho, therefore both of features and structures exhibit a ‘Crocodile’ like structure and are most probably the remnants of the amalgamation of the NQB and the NCB. The transparent reflection Moho beneath the southern part of the WG may indicate the existence of a magma channel. The Weihe Fault is interpreted as a shallow, near-surface feature resulted from the upwelling magma; SOB represents a relatively weak region and could accommodate the crustal shortening during the formation of the China continent in Triassic.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉Large-amplitude collar wave covering formation signals is still a tough problem in acoustic logging-while-drilling (LWD) measurements. In this study, we investigate the propagation and energy radiation characteristics of the monopole collar wave and the effects of grooves on reducing the interference to formation waves by finite-difference calculations. We found that the collar wave radiates significant energy into the formation by comparing the waveforms between a collar within an infinite fluid, and the acoustic LWD in different formations with either an intact or a truncated collar. The collar wave recorded on the outer surface of the collar consists of the outward-radiated energy direct from the collar (direct collar wave) and that reflected back from the borehole wall (reflected collar wave). All these indicate that the significant effects of the borehole-formation structure on collar wave were underestimated in previous studies. From the simulations of acoustic LWD with a grooved collar, we found that grooves broaden the frequency region of low collar-wave excitation and attenuate most of the energy of the interference waves by multireflections. However, grooves extend the duration of the collar wave and convert part of the collar-wave energy originally kept in the collar into long-duration Stoneley wave. Interior grooves are preferable to exterior ones because both the low-frequency and the high-frequency parts of the collar wave can be reduced and the converted inner Stoneley wave is relatively difficult to be recorded on the outer surface of the collar. Deeper grooves weaken the collar wave more greatly, but they result in larger converted Stoneley wave especially for the exterior ones. The interference waves, not only the direct collar wave but also the reflected collar wave and the converted Stoneley waves, should be overall considered for tool design.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉The oceanic crustal and uppermost lithospheric mantle structure across the Gloria Fault (GF) transcurrent plate boundary between Africa and Eurasia in the Northeast Atlantic is investigated based on seismic reflection, seismic refraction and wide-angle reflection data. This experiment used 18 ocean bottom stations along an N–S 150 km long traverse together with acquisition of a multichannel seismic reflection profile. Modeling of 〈span style="font-style:italic;"〉P〈/span〉 and 〈span style="font-style:italic;"〉S〈/span〉 seismic waves and gravimetric anomalies allowed estimation of 〈span style="font-style:italic;"〉P〈/span〉- and 〈span style="font-style:italic;"〉S〈/span〉-wave velocities, density, Poisson's ratio and discussion of a compositional model. A five-layer model is proposed in which layers 1–3 correspond to normal sediments through typical oceanic crust layers 2 and 3. Layer 5 yielded mantle velocities above 7.9 km s〈sup〉−1〈/sup〉. Layer 4 with 4 km of thickness has 〈span style="font-style:italic;"〉Vp〈/span〉 velocities between 7.1 and 7.4 km s〈sup〉−1〈/sup〉 and is clearly separated from typical oceanic crust and mantle layers. Comparison with natural analogues and published lab measurements suggest that layer 4 can be a mix of lithologies that comply with the estimated 〈span style="font-style:italic;"〉P〈/span〉 and 〈span style="font-style:italic;"〉S〈/span〉 velocities and computed Poisson's ratio and densities, such as, olivine cumulates, peridotite, gabbro and hydrated mantle. We favour the tectonic process that produces secondary porosity from which results serpentinization due to sea water circulation in fractures. Structural and seismic stratigraphic interpretation of the reflection profile shows that Neogene to recent tectonic deformation on this segment of the plate boundary concentrated on the southern side of the GF, that is, the Africa plate.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉In a lossy medium with complex frequency-dependent wave speed both rays and plane waves at an interface should satisfy the dispersion relation (that is, the wave equation), the radiation condition (the amplitude should go to zero at infinity) and the horizontal complex slowness should be continuous (Snell's law). It is known that this may lead to a transmitted wave which violates the radiation condition and which also causes problems with the phase of the reflection coefficient. In fact, ray-tracing algorithms and analytical evaluations of the reflection and transmission coefficients in anelastic media may lead to non-physical solutions related to the complex square roots of the vertical slowness and polarizations. The steepest-descent approximation with complex horizontal slowness involves non-physical complex horizontal distances, and in some cases also a non-physical vertical slowness that violates the radiation condition. Similarly, the reflection and transmission coefficients and ray-tracing codes obtained with this approach yields wrong results. In order to tackle this problem, we choose the stationary-phase approximation with real horizontal slowness. This gives real horizontal distances, the radiation condition is always satisfied and the reflection and transmission coefficients are correct. This is shown by comparison to full-wave space-time modelling results by computing the reflection and transmission coefficients and respective phase angles from synthetic seismograms. This numerical evaluation is based on a 2-D wavenumber-frequency Fourier transform. The results indicate that the stationary-phase method with a real horizontal slowness provides the correct physical solution.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉The Geysers geothermal reservoir in northern California is the site of numerous studies of both seismicity induced by injection of fluids and seismicity triggered by other earthquakes. Data from a controlled experiment in the northwest part of The Geysers in the time period 2011 to 2015 are used to study these induced and triggered earthquakes and possible differences between them. Causal solutions to the elastic equations for a porous medium show how fluid injection generates fast elastic and diffusion waves followed by a much slower diffusive wake. Calculations of fluid increment, fluid pressure and elastic stress are used to investigate both when and why seismic failure takes place. Taking into account stress concentrations caused by material heterogeneity leads to the conclusion that fluid injection by itself can cause seismic activity with no need for tectonic forces. Induced events that occur at early times are best explained by changes in stress rate, while those that occur at later times are best explained by changes in stress. While some of the seismic activity is clearly induced by injection of fluids, also present is triggered seismicity that includes aftershock sequences, swarms of seismicity triggered by other earthquakes at The Geysers and clusters of multiple earthquakes. No basic differences are found between the source mechanisms of these different types of earthquakes.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉To describe errors in the data, Gaussian distributions naturally come to mind. In many practical instances, indeed, Gaussian distributions are appropriate. In the broad field of geomagnetism, however, it has repeatedly been noted that residuals between data and models often display much sharper distributions, sometimes better described by a Laplace distribution. In this study, we make the case that such non-Gaussian behaviours are very likely the result of what is known as mixture of distributions in the statistical literature. Mixtures arise as soon as the data do not follow a common distribution or are not properly normalized, the resulting global distribution being a mix of the various distributions followed by subsets of the data, or even individual datum. We provide examples of the way such mixtures can lead to distributions that are much sharper than Gaussian distributions and discuss the reasons why such mixtures are likely the cause of the non-Gaussian distributions observed in geomagnetism. We also show that when properly selecting subdata sets based on geophysical criteria, statistical mixture can sometimes be avoided and much more Gaussian behaviours recovered. We conclude with some general recommendations and point out that although statistical mixture always tends to sharpen the resulting distribution, it does not necessarily lead to a Laplacian distribution. This needs to be taken into account when dealing with such non-Gaussian distributions.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉This paper outlines an analytical model of crack growth induced permeability changes. A theoretical solution of effective permeability of cracked porous media is derived. The fluid flow obeys Poisseuille's law along the crack and Darcy's law in the porous matrix. This solution exhibits a percolation threshold for any type of crack distribution apart from a parallel crack distribution. The physical behaviour of fluid flow through a cracked porous material is well reproduced by the proposed model. The presence of this effective permeability coupling to analytical expression of crack growth under compression enables the modelling of the permeability variation due to stress-induced cracking in a porous rock. This incorporation allows the prediction of the permeability change of a porous rock embedding an anisotropic crack distribution from any initial crack density, that is, lower, around or upper to percolation threshold. The interaction between cracks is not explicitly taken into account. The model is well applicable both to micro- and macrocracks.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Electrical conductance 〈span style="font-style:italic;"〉G〈/span〉 at 100 kHz of Berea sandstone initially saturated with varying NaCl concentrations was measured by an impedance meter at decreasing water saturation. The obtained conductance 〈span style="font-style:italic;"〉G〈/span〉 values can be well simulated by the model equation composed of conductance of bulk pore water and that of mineral surfaces by introducing both tortuosities of bulk pore water τ〈sub〉b〈/sub〉 and mineral surfaces τ〈sub〉s〈/sub〉. The surface conductivity Σ〈sub〉s〈/sub〉 = 2.1 × 10〈sup〉− 10〈/sup〉 S and the tortuosity of mineral surfaces τ〈sub〉s〈/sub〉 = 2.6 in this equation can be valid for most of the data at varying water saturation except for the lowest water saturation (〈span style="font-style:italic;"〉S〈/span〉〈sub〉w〈/sub〉 = 0.05). The tortuosity of pore water τ〈sub〉b〈/sub〉 increased from 1.7 at 〈span style="font-style:italic;"〉S〈/span〉〈sub〉w〈/sub〉 = 1.0 to 15 at 〈span style="font-style:italic;"〉S〈/span〉〈sub〉w〈/sub〉 = 0.05 with a power law relationship. The present electrical conduction model with double tortuosities of bulk pore water τ〈sub〉b〈/sub〉 and mineral surfaces τ〈sub〉s〈/sub〉 can be considered as an alternative expression of the combined Archie's first and second laws in terms of tortuosities and would be useful for describing conductance of electrolyte containing partially saturated rocks including very low water saturation.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We present the first 3 months of aftershock activity following the 2015 April 25 Gorkha earthquake 〈span style="font-style:italic;"〉M〈/span〉〈sub〉w〈/sub〉 7.8 recorded on the Nepalese Seismic network. We deployed an automatic procedure composed of three main stages: (1) coarse determination of the 〈span style="font-style:italic;"〉P〈/span〉 and 〈span style="font-style:italic;"〉S〈/span〉 onsets; (2) phase association to declare events and (3) iterative addition and refinement of onsets using the Kurtosis characteristic function. In total 9188 events could be located in the Kathmandu region with the majority having small location errors (〈4.5, 9 and 10 km in the 〈span style="font-style:italic;"〉X-, Y-〈/span〉 and 〈span style="font-style:italic;"〉Z〈/span〉-directions, respectively). Additionally, we propose a new attenuation law to estimate local magnitudes in the region. This new seismic catalogue reveals a detailed insight into the Gorkha aftershock sequence and its relation to the main shock rupture models and tectonic structures in the region. Most aftershocks fall within the Main Himalayan Thrust (MHT) shear zone or in its hangingwall. Significant temporal and lateral variations of aftershocks location are observed among them: (1) three distinct stages, highlighting subsequent jump-offs at the easternmost termination, (2) the existence of a seismic gap north of Kathmandu which matches with a low slip zone in the rupture area of the main shock, (3) the confinement of seismic activity in the trace of the May 12 〈span style="font-style:italic;"〉M〈/span〉〈sub〉w〈/sub〉 7.3 earthquake within the MHT and its hangingwall through a 30 × 30 km〈sup〉2〈/sup〉 region and (4) a shallow westward-dipping structure east of the Kathmandu klippe. These new observations with the inferred tectonic structures at depth suggest a tectonic control of part of the aftershock activity by the lateral breaks along the MHT and by the geometry of the duplex above the thrust.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉Convection in Earth's core is a viable mechanism for generating MAC waves when the top of the core is stably stratified. We quantify the generation mechanism by extending the physical description of MAC waves to include a source term due to buoyancy forces in the convecting part of the core. Solutions for the forced motion are obtained using a Green's function, which is constructed from the eigenfunctions for the unforced motion. When the source term is evaluated using the output of a numerical geodynamo model, the largest excitation occurs at even spherical harmonic degrees, corresponding to waves with symmetric azimuthal flow about the equator. We also find that the magnitude of the source term decreases at periods shorter than about 60 yr. As a result most of the wave generation is confined to waves with periods of 60 yr or longer. Quantitative predictions for the wave amplitudes depend on the projection of the source term into the eigenfunction of the waves. Strong stratification limits the penetration of density anomalies into the stratified layer, which means that the source term is confined to the lowermost part of the layer. Overtones of MAC waves with large amplitudes in the lower part of the stratified layer are more effectively generated by convection, even though these waves are heavily damped by magnetic diffusion. Generation of MAC waves by convection establishes a physical link between observable wave motion and deeper convective processes. Detection of changes in the amplitude and phase of MAC waves would constrain the generation processes and offer insights into the nature of the convection.〈/span〉

Description: 〈span class="paragraphSection"〉The paper ‘Hydrostratigraphy characterization of the Floridan aquifer system using ambient seismic noise’ was originally published in 209(2), 876–889. The original article has been updated to correct a spelling mistake in one of the authors’ names. The publisher wishes to apologise for the mistake.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The 3-D subsurface structure beneath the southern Korean Peninsula is poorly known, even though such information could be key in verifying or rejecting several competing models of the tectonic evolution of East Asia. We constructed a 3-D velocity model of the upper crust beneath the southern Korean Peninsula using 19 935 〈span style="font-style:italic;"〉P〈/span〉-wave arrivals from 747 earthquakes recorded by high-density local seismic networks. Results show significant lateral and vertical variations: velocity increases from northwest to southeast at shallow depths, and significant velocity variations are observed across the South Korea Tectonic Line between the Okcheon Fold Belt and the Youngnam Massif. Collision between the North and South China blocks during the Early Cretaceous might have caused extensive deformation and the observed negative velocity anomalies in the region. The results of the tomographic inversion, combined with the findings of previous studies of Bouguer and isostatic gravity anomalies, indicate the presence of high-density material in the upper and middle crust beneath the Gyeongsang Basin in the southeastern Korean Peninsula. Although our results partially support the indentation tectonic model, it is still premature to discard other tectonic evolution models because our study only covers the southern half of the peninsula.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The perfectly matched layer (PML) is an efficient absorbing technique for numerical wave simulation. The complex frequency-shifted PML (CFS-PML) introduces two additional parameters in the stretching function to make the absorption frequency dependent. This can help to suppress converted evanescent waves from near grazing incident waves, but does not efficiently absorb low-frequency waves below the cut-off frequency. To absorb both the evanescent wave and the low-frequency wave, the double-pole CFS-PML having two poles in the coordinate stretching function was developed in computational electromagnetism. Several studies have investigated the performance of the double-pole CFS-PML for seismic wave simulations in the case of a narrowband seismic wavelet and did not find significant difference comparing to the CFS-PML. Another difficulty to apply the double-pole CFS-PML for real problems is that a practical strategy to set optimal parameter values has not been established. In this work, we study the performance of the double-pole CFS-PML for broad-band seismic wave simulation. We find that when the maximum to minimum frequency ratio is larger than 16, the CFS-PML will either fail to suppress the converted evanescent waves for grazing incident waves, or produce visible low-frequency reflection, depending on the value of α. In contrast, the double-pole CFS-PML can simultaneously suppress the converted evanescent waves and avoid low-frequency reflections with proper parameter values. We analyse the different roles of the double-pole CFS-PML parameters and propose optimal selections of these parameters. Numerical tests show that the double-pole CFS-PML with the optimal parameters can generate satisfactory results for broad-band seismic wave simulations.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Density variations drive mass transport in the Earth from plate tectonics to convection in the mantle and core. Nevertheless, density remains poorly known because most geophysical measurements used to probe the Earth's interior either have little sensitivity to density, suffer from trade-offs or from non-uniqueness. With the ongoing expansion of computational power, it has become possible to accurately model complete seismic wavefields in a 3-D heterogeneous Earth, and to develop waveform inversion techniques that account for complicated wavefield effects. This may help to improve resolution of density. Here, we present a pilot study where we explore the extent to which waveform inversion may be used to better recover density as a separate, independent parameter. We perform numerical simulations in 2-D to investigate under which conditions, and to what extent density anomalies may be recovered in the Earth's mantle. We conclude that density can indeed be constrained by seismic waveforms, mainly as a result of scattering effects at density contrasts. As a consequence, the low-frequency part of the wavefield is the most important for constraining the actual extent of anomalies. While the impact of density heterogeneities on the wavefield is small compared to the effects of velocity variations, it is likely to be detectable in modern regional- to global-scale measurements. We also conclude that the use of gravity data as additional information does not help to further improve the recovery of density anomalies unless strong 〈span style="font-style:italic;"〉a priori〈/span〉 constraints on the geometry of density variations are applied. This is a result of the inherent physical non-uniqueness of potential-field inverse problems. Finally, in the limited numerical setup that we employ, we find that the initially supplied anomalies in 〈span style="font-style:italic;"〉S〈/span〉- and 〈span style="font-style:italic;"〉P〈/span〉-velocity models are of minor importance.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The accuracy of seismic numerical simulations and the effectiveness of imaging conditions are important in reverse time migration studies. Using the pseudospectral method, the precision of the calculated spatial derivative of the seismic wavefield can be improved, increasing the vertical resolution of images. Low-frequency background noise, generated by the zero-lag cross-correlation of mismatched forward-propagated and backward-propagated wavefields at the impedance interfaces, can be eliminated effectively by using the imaging condition based on the wavefield decomposition technique. The computation complexity can be reduced when imaging is performed in the frequency domain. Since the Fourier transformation in the 〈span style="font-style:italic;"〉z〈/span〉-axis may be derived directly as one of the intermediate results of the spatial derivative calculation, the computation load of the wavefield decomposition can be reduced, improving the computation efficiency of imaging. Comparison of the results for a pulse response in a constant-velocity medium indicates that, compared with the finite difference method, the peak frequency of the Ricker wavelet can be increased by 10–15 Hz for avoiding spatial numerical dispersion, when the second-order spatial derivative of the seismic wavefield is obtained using the pseudospectral method. The results for the SEG/EAGE and Sigsbee2b models show that the signal-to-noise ratio of the profile and the imaging quality of the boundaries of the salt dome migrated using the pseudospectral method are better than those obtained using the finite difference method.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We present an algorithm for efficient 3-D inversion of marine controlled-source electromagnetic data. The efficiency is achieved by exploiting the redundancy in data. The data redundancy is reduced by compressing the data through stacking of the response of transmitters which are in close proximity. This stacking is equivalent to synthesizing the data as if the multiple transmitters are simultaneously active. The redundancy in data, arising due to close transmitter spacing, has been studied through singular value analysis of the Jacobian formed in 1-D inversion. This study reveals that the transmitter spacing of 100 m, typically used in marine data acquisition, does result in redundancy in the data. In the proposed algorithm, the data are compressed through stacking which leads to both computational advantage and reduction in noise. The performance of the algorithm for noisy data is demonstrated through the studies on two types of noise, viz., uncorrelated additive noise and correlated non-additive noise. It is observed that in case of uncorrelated additive noise, up to a moderately high (10 percent) noise level the algorithm addresses the noise as effectively as the traditional full data inversion. However, when the noise level in the data is high (20 percent), the algorithm outperforms the traditional full data inversion in terms of data misfit. Similar results are obtained in case of correlated non-additive noise and the algorithm performs better if the level of noise is high. The inversion results of a real field data set are also presented to demonstrate the robustness of the algorithm. The significant computational advantage in all cases presented makes this algorithm a better choice.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The Kaapvaal craton (South Africa) was the host of several major magmatic events during the Palaeoproterozoic, including the volcanic Hekpoort and Ongeluk Formations. Their possible comagmatic origin is the subject of a long debate. We performed a palaeomagnetic study of the Hekpoort Formation to be compared with the available palaeopole of the Ongeluk Formation, but also to contribute to the apparent polar wander path of the Kaapvaal craton. Characterization of magnetic mineralogy by three-axis thermal demagnetization of isothermal remanent magnetization and magnetic susceptibility versus temperature points out magnetite as the main remanence carrier in most samples.Five magnetic components were identified in total, of which the least stable (HKE) near parallels the present geomagnetic field. At higher levels of demagnetization (above 400 °C), two components (HKD and HKC) are identified as thermoviscous overprints likely related to the Karoo large igneous province (LIP) and a magmatic event which occurred between the emplacement of the ∼2055 Ma Bushveld Complex and HKD (possibly linked to the Umkondo LIP), respectively. This LIP is known to be associated with extensive remagnetization. The second most stable component HKB was also revealed at higher steps of thermal and alternative-field treatment. The HKB palaeopole (latitude = 28.4°N and longitude = 54°E) is similar to those reported from the Bushveld Complex (∼2055 Ma) and the Vredefort impact structure (∼2023 Ma). A potentially primary remanence direction (HKA; declination = 337°, inclination = 80° and α = 6.2°) was identified in most sites during the highest levels of thermal demagnetization. Note that the HKA pole position (latitude = −44°N and longitude = 40°E) is significantly different from the palaeopole for the Ongeluk Formation (latitude = −0.5°N and longitude = 107°E). Although, the primary nature of HKA is supported by positive fold and reversal tests, we cannot exclude the possibility that this component represents an overprint. HKA is, however, most likely older than ∼2.0 Ga given its anteriority to HKB components.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Displacement spectra from accelerograms recorded within a few kilometres (〈3.5 km) of nuclear and chemical explosion shotpoints indicate that the amplitude spectral level is not flat below the source corner frequency (〈span style="font-style:italic;"〉fc〈/span〉), but increases gradually towards zero frequency. During an explosion, the volume in the immediate vicinity of the shotpoint is inelastic. In this paper, we develop a time-domain expression for the deformation at the elastic boundary limit of this volume, which in the frequency domain supports this observation. We refer to this expression as the time-domain source function (TDSF) of an explosion. The proposed TDSF has two terms, the first term representing a ‘static’ contribution and the second term a ‘dynamic’ contribution to the total deformation field. The static contribution dominates over the dynamic contribution at frequencies below 〈span style="font-style:italic;"〉fc〈/span〉 and causes the gradual increase in the spectral level. For the low-yield under/or overburied explosions (yield 〈5 Kt), 〈span style="font-style:italic;"〉fc〈/span〉 is relatively high and this increase in spectral amplitude is pronouncedly observed. The correct interpretation of the observed spectral amplitudes below 〈span style="font-style:italic;"〉fc〈/span〉 can, therefore, play a crucial role in estimating source parameters of explosions. For 〈span style="font-style:italic;"〉f〈/span〉 〉 〈span style="font-style:italic;"〉fc〈/span〉, the dynamic contributions dominate and decay approximately as 〈span style="font-style:italic;"〉f〈/span〉〈sup〉 − 2〈/sup〉. For seismic waves propagating from the boundary 〈span style="font-style:italic;"〉R〈/span〉〈sub〉el〈/sub〉 (a transition limit of the non-linear to the elastic zone) to large distances, the static and the dynamic wavefields are affected identically by attenuation and spreading. Hence, the attenuation corrected large distance explosion spectra should exhibit these spectral characteristics. An analysis of the regional 〈span style="font-style:italic;"〉P〈/span〉-wave spectra from a few low-yield explosions provides evidence for this finding and also for the yield scaling by a factor of 2 between the nuclear and chemical explosions, especially for similar emplacement conditions. We also illustrate that when convolved with a time function [exp( − 〈span style="font-style:italic;"〉C〈/span〉/〈span style="font-style:italic;"〉R〈/span〉〈sub〉el〈/sub〉) 〈span style="font-style:italic;"〉H〈/span〉(〈span style="font-style:italic;"〉t〈/span〉)], where 〈span style="font-style:italic;"〉C〈/span〉 is the material velocity at the shotpoint, the TDSF yields the reduced displacement potential (RDP) at 〈span style="font-style:italic;"〉R〈/span〉〈sub〉el〈/sub〉 of the source. We use this proposed RDP to investigate the influence of yield and depth of burial on the spectral overshoot and 〈span style="font-style:italic;"〉fc〈/span〉 of explosion sources.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉After more than 13 yr of GRACE monthly data, the determined secular trend of gravity field variation can be used to study the regions of glacial isostatic adjustment (GIA). Here we focus on Fennoscandia where long-term terrestrial and high-quality GPS data are available, and we study the monthly GRACE data from three analysis centres. We present a new approximate formula to convert the secular trend of the GRACE gravity change to the land uplift rate without making assumptions of the ice load history. The question is whether the GRACE-derived land uplift rate by our method is related to GIA. A suitable post-processing method for the GRACE data is selected based on weighted RMS differences with the GPS data. The study reveals that none of the assumed periodic changes of the GRACE gravity field is significant in the estimation of the secular trend, and they can, therefore, be neglected. Finally, the GRACE-derived land uplift rates are obtained using the selected post-processing method, and they are compared with GPS land uplift rate data. The GPS stations with significant differences were marked using a statistical significance test. The smallest rms difference (1.0 mm a〈sup〉–1〈/sup〉) was obtained by using GRACE data from the University of Texas.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Seismic monitoring provides valuable information regarding the time-varying changes in subsurface physical properties caused by natural or man-made processes. However, the resulting changes in the earth's subsurface properties are often small both in terms of magnitude and spatial extent, leading to minimal time-lapse differences in seismic amplitude or traveltime. In order to better extract information from the time-lapse data, we show that exploiting the full seismic waveform information can be critical. In this study, we develop and test methods of full waveform inversion that estimate an optimal subsurface model of time-varying elastic properties in order to fit the observed time-lapse seismic data with predicted waveforms based on numerical solutions of the wave equation. Time-lapse full waveform inversion is nonlinear and non-unique, and depends on the knowledge of the baseline velocity model before a change, and (non-)repeatability of earthquake source and sensor parameters, and of ambient and cultural noise. We propose to use repeating earthquake data sets acquired with permanent arrays of seismic sensors to enhance the repeatability of source and sensor parameters. We further develop and test time-lapse parallel, double-difference and bootstrapping inversion strategies to mitigate the dependence on the baseline velocity model. The parallel approach uses a time-invariant full waveform inversion method to estimate velocity models independently of the different source event times. The double-difference approach directly estimates velocity changes from time-lapse waveform differences, requiring excellent repeatability. The bootstrapping approach inverts for velocity models sequentially in time, implicitly constraining the time-lapse inversions, while relaxing an explicit requirement for high data repeatability. We assume that prior to the time-lapse inversion, we can estimate the true source locations and the origin time of the events, and also we can also obtain a reasonably accurate baseline velocity model. Extensive synthetic tests using a realistic velocity model developed from a real project area demonstrate the potential of full waveform inversion to estimate velocity changes from dense surface arrays of seismic stations recording a small number of repeating events. Analysis of sensitivity kernels suggests that positioning sensors at large distances allows for a stable recovery of velocity changes near the event locations by illuminating the inversion area with a wide aperture of angles. We show how full waveform inversion maps the errors in the baseline velocity model and the non-repeatable noise into the estimates of time-lapse velocity changes. Among the three time-lapse inversion methods, parallel inversion is most affected by non-repeatability factors, and is thus the least robust and most contaminated by artefacts. In contrast, the double-difference and bootstrapping methods result in more accurate time-lapse inversions. As the non-repeatability of both sources and noise increases, the bootstrapping method provides more robust and accurate results than the double-difference method.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We investigated a new technique for aquifer characterization that uses cross-correlation of ambient seismic noise to determine seismic velocity structure of the Floridan aquifer system (FAS). Accurate characterization of aquifer systems is vital to hydrogeological research and groundwater management but is difficult due to limited subsurface data and heterogeneity. Previous research on the carbonate FAS found that confining units and high permeability flow zones have distinct seismic velocities. We deployed an array of 9 short period seismometers from 11/2013 to 3/2014 in Indian Lake State Forest near Ocala, Florida, to image the hydrostratigraphy of the aquifer system using ambient seismic noise. We find that interstation distance strongly influences the upper and lower frequency limits of the data set. Seismic waves propagating within 1.5 and 7 wavelengths between stations were optimal for reliable group velocity measurements and both an upper and lower wavelength threshold was used. A minimum of 100–250 hr of signal was needed to maximize signal-to-noise ratio and to allow cross-correlation convergence. We averaged measurements of group velocity between station pairs at each frequency band to create a network average dispersion curve. A family of 1-D shear-wave velocity profiles that best represents the network average dispersion was then generated using a Markov Chain Monte Carlo (MCMC) algorithm. The MCMC algorithm was implemented with either a fixed number of layers, or as transdimensional in which the number of layers was a free parameter. Results from both algorithms require a prominent velocity increase at ∼200 m depth. A shallower velocity increase at ∼60 m depth was also observed, but only in model ensembles created by collecting models with the lowest overall misfit to the observed data. A final round of modelling with additional prior constraints based on initial results and well logs produced a mean shear-wave velocity profile taken as the preferred solution for the study site. The velocity increases at ∼200 and ∼60 m depth are consistent with the top surfaces of two semi-confining units of the study area and the depths of high-resistivity dolomite units seen in geophysical logs and cores from the study site. Our results suggest that correlation of ambient seismic noise holds promise for hydrogeological investigations. However, complexities in the cross-correlations at high frequencies and short traveltimes at low frequencies added uncertainty to the data set.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉This paper is concerned with the applicability of Pareto Multi-Objective Global Optimization (PMOGO) algorithms for solving different types of geophysical inverse problems. The standard deterministic approach is to combine the multiple objective functions (i.e. data misfit, regularization and joint coupling terms) in a weighted-sum aggregate objective function and minimize using local (decent-based) smooth optimization methods. This approach has some disadvantages: (1) appropriate weights must be determined for the aggregate, (2) the objective functions must be differentiable and (3) local minima entrapment may occur. PMOGO algorithms can overcome these drawbacks but introduce increased computational effort. Previous work has demonstrated how PMOGO algorithms can overcome the first issue for single data set geophysical inversion, that is, the trade-off between data misfit and model regularization. However, joint inversion, which can involve many weights in the aggregate, has seen little study. The advantage of PMOGO algorithms for the other two issues has yet to be addressed in the context of geophysical inversion. In this paper, we implement a PMOGO genetic algorithm and apply it to physical-property-, lithology- and surface-geometry-based inverse problems to demonstrate the advantages of using a global optimization strategy. Lithological inversions work on a mesh but use integer model parameters representing rock unit identifiers instead of continuous physical properties. Surface geometry inversions change the geometry of wireframe surfaces that represent the contacts between discrete rock units. Despite the potentially high computational requirements of global optimization algorithms (compared to local), their application to realistically sized 2-D geophysical inverse problems is within reach of current capacity of standard computers. Furthermore, they open the door to geophysical inverse problems that could not otherwise be considered through traditional optimization strategies.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉We present a skeletonized inversion method that inverts surface-wave data for the 〈span style="font-style:italic;"〉Qs〈/span〉 quality factor. Similar to the inversion of dispersion curves for the 〈span style="font-style:italic;"〉S〈/span〉-wave velocity model, the complicated surface-wave arrivals are skeletonized as simpler data, namely the amplitude spectra of the windowed Rayleigh-wave arrivals. The optimal 〈span style="font-style:italic;"〉Qs〈/span〉 model is the one that minimizes the difference in the peak frequencies of the predicted and observed Rayleigh wave arrivals using a gradient-based wave-equation optimization method. Solutions to the viscoelastic wave-equation are used to compute the predicted Rayleigh-wave arrivals and the misfit gradient at every iteration. This procedure, denoted as wave-equation 〈span style="font-style:italic;"〉Qs〈/span〉 inversion (WQ〈sub〉〈span style="font-style:italic;"〉s〈/span〉〈/sub〉), does not require the assumption of a layered model and tends to have fast and robust convergence compared to full waveform inversion (FWI). Numerical examples with synthetic and field data demonstrate that the WQ〈sub〉〈span style="font-style:italic;"〉s〈/span〉〈/sub〉 method can accurately invert for a smoothed approximation to the subsurface 〈span style="font-style:italic;"〉Qs〈/span〉 distribution as long as the 〈span style="font-style:italic;"〉Vs〈/span〉 model is known with sufficient accuracy.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The International Association of Geodesy released in July 2015 a resolution for the definition and realisation of an International Height Reference System (IHRS). According to this resolution, the IHRS coordinates are potential differences referring to the equipotential surface of the Earth's gravity field realised by the conventional value 〈span style="font-style:italic;"〉W〈/span〉〈sub〉0〈/sub〉 = 62 636 853.4 m〈sup〉2〈/sup〉s〈sup〉−2〈/sup〉. A main component of the IHRS realisation is the integration of the existing height systems into the global one; that is existing vertical coordinates should be referred to one and the same reference level realised by the conventional 〈span style="font-style:italic;"〉W〈/span〉〈sub〉0〈/sub〉. This procedure is known as vertical datum unification and its main result are the vertical datum parameters, that is the potential differences between the local and the global reference levels. In this paper, we rigorously derive the observation equations for the vertical datum unification in terms of potential quantities based on the geodetic boundary value problem (GBVP) approach. Those observation equations are then empirically evaluated for the vertical datum unification of the North American and South American height systems. In the first case, simulations performed in North America provide numerical estimates about the impact of omission errors and direct and indirect effects on the vertical datum parameters. In the second case, a combination of local geopotential numbers, ITRF coordinates, satellite altimetry observations, tide gauge registrations and high-resolution gravity field models is performed to estimate the level differences between the South American height systems and the global level 〈span style="font-style:italic;"〉W〈/span〉〈sub〉0〈/sub〉. Results show that indirect effects vanish when a satellite-only gravity field model with a degree higher than 〈span style="font-style:italic;"〉n〈/span〉 ≥ 180 is used for the solution of the GBVP. However, the component derived from satellite-only global gravity models has to be refined with terrestrial gravity data to minimise the omission error and its effect on the vertical datum parameter estimation. The empirical evaluations demonstrate that the vertical datum unification should be based on geodetic stations of highest quality and standardised geodetic data; for example, geometric coordinates should refer to the same ITRF and be given in the same tide system and reference epoch as the geopotential numbers and gravity field model. After a standardisation of the input data used in the unification of the South American height systems and a rigorous error propagation analysis, we demonstrate that the vertical datum parameters can be estimated with accuracy better than ±5 cm in well-surveyed regions and some decimetres (±40 cm) in sparsely surveyed regions. This paper concludes with detailed guidelines for the appropriate data treatment when the integration of a local vertical datum into the IHRS is desired. These guidelines may be applicable in any region of the world.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉In order to upscale the induced polarization (IP) response of porous media, from the pore scale to the sample scale, we implement a procedure to compute the macroscopic complex resistivity response of random tube networks. A network is made of a 2-D square-meshed grid of connected tubes, which obey to a given tube radius distribution. In a simplified approach, the electrical impedance of each tube follows a local Pelton resistivity model, with identical resistivity, chargeability and Cole–Cole exponent values for all the tubes—only the time constant varies, as it depends on the radius of each tube and on a diffusion coefficient also identical for all the tubes. By solving the conservation law for the electrical charge, the macroscopic IP response of the network is obtained. We fit successfully the macroscopic complex resistivity also by a Pelton resistivity model. Simulations on uncorrelated and correlated networks, for which the tube radius distribution is so that the decimal logarithm of the radius is normally distributed, evidence that the local and macroscopic model parameters are the same, except the Cole–Cole exponent: its macroscopic value diminishes with increasing heterogeneity (i.e. with increasing standard deviation of the radius distribution), compared to its local value. The methodology is also applied to six siliciclastic rock samples, for which the pore radius distributions from mercury porosimetry are available. These samples exhibit the same behaviour as synthetic media, that is, the macroscopic Cole–Cole exponent is always lower than the local one. As a conclusion, the pore network method seems to be a promising tool for studying the upscaling of the IP response of porous media.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The extent to which aseismic deformation relaxes co-seismic stress changes on a fault zone is fundamental to assessing the future seismic hazard following any earthquake, and in understanding the mechanical behaviour of faults. Here we use models of stress-driven afterslip and viscoelastic relaxation, in conjunction with post-seismic InSAR measurements, to show that there has been minimal release of co-seismic stress changes through post-seismic deformation following the 2003 〈span style="font-style:italic;"〉M〈/span〉〈sub〉w〈/sub〉 6.6 Bam earthquake. Our analysis indicates the faults at Bam remain predominantly locked, suggesting that the co- plus interseismically accumulated elastic strain stored downdip of the 2003 rupture patch may be released in a future 〈span style="font-style:italic;"〉M〈/span〉〈sub〉w〈/sub〉 6 earthquake. Our observations and models also provide an opportunity to probe the growth of topography at Bam. We find that, for our modelled afterslip distribution to be consistent with forming the sharp step in the local topography over repeated earthquake cycles, and also to be consistent with the geodetic observations, requires either (1) far-field tectonic loading equivalent to a 2–10 MPa deviatoric stress acting across the fault system, which suggests it supports stresses 60–100 times less than classical views of static fault strength, or (2) that the fault surface has some form of mechanical anisotropy, potentially related to corrugations on the fault plane, that controls the sense of slip.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉The Deep Fault Drilling Project (DFDP) on the central Alpine Fault, South Island, New Zealand, has motivated a broad range of geophysical and geological studies intended to characterize the fault system in the locality of the drill site at various scales. In order to better understand the structural features of the central Alpine Fault, we have developed 3-D 〈span style="font-style:italic;"〉P〈/span〉- and 〈span style="font-style:italic;"〉S〈/span〉-wave velocity (〈span style="font-style:italic;"〉V〈/span〉〈sub〉P〈/sub〉 and 〈span style="font-style:italic;"〉V〈/span〉〈sub〉S〈/sub〉) models of the region by double-difference tomography using data sets from multiple seismic networks. In previous work, the quality of the 〈span style="font-style:italic;"〉S〈/span〉-wave model has been poor due to the small number of available 〈span style="font-style:italic;"〉S〈/span〉-wave picks. We have utilized a new high-accuracy automatic 〈span style="font-style:italic;"〉S〈/span〉-wave picker to increase the number of usable 〈span style="font-style:italic;"〉S〈/span〉-wave arrivals by more than a factor of two, thereby substantially improving the 〈span style="font-style:italic;"〉V〈/span〉〈sub〉S〈/sub〉 model. Compared to previous studies, our new higher-resolution 〈span style="font-style:italic;"〉V〈/span〉〈sub〉P〈/sub〉 model based on more observations shows a clear 〈span style="font-style:italic;"〉V〈/span〉〈sub〉P〈/sub〉 contrast (higher 〈span style="font-style:italic;"〉V〈/span〉〈sub〉P〈/sub〉 on the southeast hanging wall side) at depths of 5–10 km near the DFDP drill sites. With our better resolved 〈span style="font-style:italic;"〉V〈/span〉〈sub〉S〈/sub〉 model, in the same region, we detect a sharply defined high 〈span style="font-style:italic;"〉V〈/span〉〈sub〉S〈/sub〉 body (〈span style="font-style:italic;"〉V〈/span〉〈sub〉S〈/sub〉 〉 3.7 km s〈sup〉−1〈/sup〉) within the hanging wall. Our earthquake relocations reveal the presence of clusters within and around low-velocity zones in the hanging wall southeast of the Alpine Fault. Together with the improved earthquake locations, the 〈span style="font-style:italic;"〉P〈/span〉- and 〈span style="font-style:italic;"〉S〈/span〉-wave tomography results reveal the Alpine Fault to be marked by a velocity contrast throughout most of the study region. The fault dips southeastwards at about 50° from 5 to 15 km depth, as inferred from the velocity structure, seismicity and observations of fault zone guided waves.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Vp/Vs ratio where Vp and Vs are 〈span style="font-style:italic;"〉P〈/span〉- and 〈span style="font-style:italic;"〉S〈/span〉-wave velocities is an indicator of rock composition, but estimates of Vp/Vs for the lower continental crust remain sparse. We present estimates of Vs, Vp and Vp/Vs in the crust of the Archean-Paleoproterozoic Siberian craton that are obtained by simultaneous inversion of 〈span style="font-style:italic;"〉P〈/span〉 and 〈span style="font-style:italic;"〉S〈/span〉 receiver functions from GSN seismograph stations NRIL, YAK and TIXI. These stations are located in the region of the Siberian traps (NRIL), close to the Laptev Sea Rift (TIXI) and the Viluy rift system (YAK). The most conspicuous result of our analysis is a high Vp/Vs ratio (≥2.0) at depths from 20–30 to 40 km. A very high Vp in this layer (from 7 to 8 km s〈sup〉−1〈/sup〉) is indicative of magmatic underplating. We find broadly similar data in the western Mediterranean and in India. In a dry lower crust the Vp/Vs ratio is ∼1.8, which is hard to reconcile with the estimates 〉2.0. A coincidence in depths of zones of high electric conductivity and of anomalously high Vp/Vs in Siberia suggests that both may have the same origin: fluid-filled porosity. The porosity which is required by our seismic observations is on the order of 1 per cent. Origins of the fluids may be linked with processes of magmatic underplating.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Abstract〈/div〉Discrepancies between geological, seismic and geodetic rates of strain can indicate that rates of crustal deformation, and hence seismic hazard, are varying through time. Previous studies in the northern Shanxi Grabens, at the northeastern corner of the Ordos Plateau in northern China, have found extension rates of anywhere between 0 and 6 mm a〈sup〉−1〈/sup〉 at an azimuth of between 95° and 180°. In this paper we determine extension rates across the northern Shanxi Grabens from offset geomorphological features and a variety of Quaternary dating techniques (including new IRSL and Ar-Ar ages), a Kostrov summation using a 700 yr catalogue of historical earthquakes, and recent campaign GPS measurements. We observe good agreement between Quaternary, seismic and geodetic rates of strain, and we find that the northern Shanxi Grabens are extending at around 1–2 mm a〈sup〉−1〈/sup〉 at an azimuth of ≈151°. The azimuth of extension is particularly well constrained and can be reliably inferred from catalogues of small earthquakes. We do not find evidence for any substantial variations in extension rate through time, though there is a notable seismic moment rate deficit since 1750. This deficit could indicate complex fault interactions across large regions, aseismic accommodation of deformation, or that we are quite late in the earthquake cycle with the potential for larger earthquakes in the relatively near future.〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉We investigate a numerical performance of four different schemes applied to a regional recovery of the gravity anomalies from the third-order gravitational tensor components (assumed to be observable in the future) synthetized at the satellite altitude of 200 km above the mean sphere. The first approach is based on applying a regional inversion without modelling the far-zone contribution or long-wavelength support. In the second approach we separate integral formulas into two parts, that is, the effects of the third-order disturbing tensor data within near and far zones. Whereas the far-zone contribution is evaluated by using existing global geopotential model (GGM) with spectral weights given by truncation error coefficients, the near-zone contribution is solved by applying a regional inversion. We then extend this approach for a smoothing procedure, in which we remove the gravitational contributions of the topographic-isostatic and atmospheric masses. Finally, we apply the remove-compute-restore (r-c-r) scheme in order to reduce the far-zone contribution by subtracting the reference (long-wavelength) gravity field, which is computed for maximum degree 80. We apply these four numerical schemes to a regional recovery of the gravity anomalies from individual components of the third-order gravitational tensor as well as from their combinations, while applying two different levels of a white noise. We validated our results with respect to gravity anomalies evaluated at the mean sphere from EGM2008 up to the degree 250. Not surprisingly, better fit in terms of standard deviation (STD) was attained using lower level of noise. The worst results were gained applying classical approach, STD values of our solution from 〈span style="font-style:italic;"〉Tzzz〈/span〉 are 1.705 mGal (noise value with a standard deviation 0.01 × 10〈sup〉 − 15〈/sup〉m〈sup〉 − 1〈/sup〉s〈sup〉 − 2〈/sup〉) and 2.005 mGal (noise value with a standard deviation 0.05 × 10〈sup〉 − 15〈/sup〉m〈sup〉 − 1〈/sup〉s〈sup〉 − 2〈/sup〉), while the superior from r-c-r up to the degree 80, STD fit of gravity anomalies from 〈span style="font-style:italic;"〉Tzzz〈/span〉 with respect to the same counterpart from EGM2008 is 0.510 mGal (noise value with a standard deviation 0.01 × 10〈sup〉 − 15〈/sup〉m〈sup〉 − 1〈/sup〉s〈sup〉 − 2〈/sup〉) and 1.190 mGal (noise value with a standard deviation 0.05 × 10〈sup〉 − 15〈/sup〉m〈sup〉 − 1〈/sup〉s〈sup〉 − 2〈/sup〉).〈/span〉

Description: 〈span class="paragraphSection"〉〈div class="boxTitle"〉Summary〈/div〉Seismic arrays provide an important means of enhancing seismic signals and determining the directional properties of the wavefield by beamforming. When multiple arrays are to be used together, the viewpoint needs to be modified from looking outwards from each array to focusing on a specific target area and so constraining the portions of the waveforms to be analysed. Beamforming for each array is supplemented by the relative time constraints for propagation from the target to each array to provide tight spatial control. Simultaneous multiple array analysis provides a powerful tool for source characterization, and for structural analysis of scatterers as virtual sources. The multiple array concept allows us to illuminate a specific point in the Earth from many different directions and thus maps detailed patterns of heterogeneity in the Earth. Furthermore, illumination of the structure from multiple directions using data from the same event minimizes source effects to provide clearer images of heterogeneity. The analysis is based on a similar concept to the backprojection technique, where a part of the seismic wave train is mapped to a specific point in space by ray tracing. In contrast to classic backprojection where the incoming energy is mapped onto a horizontal plane with limited vertical resolution, the multiarray method controls depth response by combining relative time constraints between the arrays and conventional beamforming. We illustrate this approach with application to two earthquakes at moderate depth. The results show that the use of simultaneous multiple arrays can provide improvement both in signal quality and resolution, with the additional benefit of being able to accurately locate the source of the incoming energy and map large areas with only a limited number of such arrays.〈/span〉

Description: By introducing two types of zenith troposphere delay (ZTD) products in precise point positioning (PPP), we developed the ZTD-corrected PPP and the ZTD-constrained PPP, both of them reduced the PPP convergence time. Both enhanced PPP methods are examined by global empirical ZTD models and regional ZTD corrections. For global ZTD models, we verified that ZTD-corrected PPP will deviate the positioning results, while ZTD-constrained PPP could produce unbiased estimations. Therefore, the latter is utilized to study the performance of global ZTD models (ITG, GPT2w, GZTD and UNB3m). After numerous experiments, we found that the performance of ZTD models was positively related to the real ZTD accuracy, and we proposed a universal tropospheric stochastic model 2SQR(9rms) which denotes double the square of nine times ZTD rms, to constrain ZTD in PPP. The proposed model subsequently was validated by real-time static and kinematic ZTD-constrained PPP on the premise that the ZTD rms on every station was known. Compared with traditional PPP, in static PPP, the number of improved stations is increased by 15.5 per cent (ITG), 14.4 per cent (GPT2w), 11.1 per cent (GZTD) and 8.3 per cent (UNB3m). For kinematic PPP, PPP constrained by ITG model still had the best performance, the number of improved stations is increased by 14.4 per cent, after 30 min of initialization time, 13.4 cm east, 13.4 cm north and 11.7 cm up positioning accuracy was obtained, compared with 15.3 cm east, 15.3 cm north and 14.3 cm up accuracy by traditional PPP. In addition, experiments using regional ZTD corrections to enhance real-time PPP showed that both ZTD-corrected PPP and ZTD-constrained PPP can notably reduce the convergence time on the vertical component (within 15 cm).

Description: We investigated post-seismic velocity changes within the fault zone of the 2008 M 7.9 Wenchuan earthquake using coda wave data of repeating small earthquakes. We employed template matching and grid search methods to identify well-defined repeating earthquakes in order to minimize artefacts induced by variations in source location. We identified a total of 12 isolated patches in the fault zone that ruptured more than twice in a 1 yr period after the M 7.9 earthquake. We applied the coda wave interferometry technique to the waveform data of the 34 identified repeating earthquakes to estimate velocity changes between the first and subsequent events in each cluster. We found that major post-seismic velocity changes occurred in the southwestern part of the rupture area, where the main rupture was initiated and characterized by thrust motion, while the Beichuan area in the northeastern part of the rupture zone appears to experience very little post-seismic velocity changes.

Description: In this study, we demonstrate the feasibility of imaging broad-band (10–150 s) Rayleigh wave phase velocity maps on a continental scale using ambient noise tomography (ANT). We obtain broad-band Rayleigh waves from cross-correlations of ambient noise data between all station pairs of USArray and measure the dispersion curves from these cross-correlations at a period band of 10–150 s. The large-scale dense USArray enables us to obtain over 500 000 surface wave paths which cover the contiguous United States densely. Using these paths, we generate Rayleigh wave phase velocity maps at 10–150 s periods. Our phase velocity maps are similar to other reported phase velocity maps based on ambient noise data at short periods (〈50 s) and based on earthquake data at intermediate/long periods (50–90 s). This study extends ANT from short/intermediate periods (〈50 s) to long periods up to 150 s in a continental scale of the USA. These broad-band phase velocity maps from ANT can be used to construct 3-D lithospheric and asthenospheric velocity structures.

Description: Clinopyroxene (Cpx) phenocrysts have the potential to record magmatic water contents and magma ascent rates through their concentration of H 2 O (incorporated as H + and commonly referred to as water). Here we investigate three issues related to the fidelity and utility of the Cpx water record: the partitioning of water between Cpx and melt, the diffusivity of water in natural Cpx phenocrysts, and the possibility for water loss in Cpx erupted in lava flows. Samples studied are from volcanic ash of the 1974 eruption of Volcán de Fuego (Guatemala) and scoria and lava from the 1977 eruption on Seguam Island (Alaska). The partitioning of water was determined by analyzing melt inclusions (MIs) and the adjacent clinopyroxene host by ion microprobe. For seven Cpx-hosted MIs from Seguam, the partition coefficient is well predicted (within 0·1 wt % H 2 O on average) using the temperature-dependent parameterization of O’Leary et al. (2010; Earth and Planetary Science Letters 297 , 111–120). For the determination of diffusivity, H 2 O concentration profiles were measured in oriented Cpx from Fuego by ion microprobe. Water decreases toward the rim, consistent with diffusive re-equilibration during ascent-driven degassing. Using previously estimated durations of ascent (7–12 min), we determined the H + diffusivity (10 –9·7 –10 –10·3 m 2 s –1 at a temperature of 1030 °C) that would satisfy these timescales. These calculated DH+Cpx values are comparable with the medial values determined for natural Cpx (Mg# 〈 92·5) in laboratory diffusion studies. Tephra and lava co-erupted on Seguam bear similar Cpx populations in their major and trace element compositions, but the lava Cpx contain 80% lower H 2 O concentrations on average than the tephra Cpx. Using the DH+Cpx values obtained from the Fuego Cpx, the difference in H 2 O between the lava and tephra Cpx can be attributed to post-eruption H 2 O loss during the estimated ~20 min the sample remained above the H + closure temperature. The results from this study indicate that clinopyroxene from slowly cooled basaltic lavas should not be used to reconstruct initial magmatic water contents. High ascent rates, rapid post-eruptive cooling, large phenocrysts, and cooler magma (e.g. andesites and rhyolites) favor better preservation of water in Cpx. In cases of H + loss, Cpx zonation can be exploited as a chronometer for syn- and post-eruptive volcanic processes.

Description: The genesis of primitive arc magmas has had a major impact on continent formation through time, but the rarity of exposures of deep arc sections limits our understanding of the details of melt migration and differentiation. Abundant pyroxenites are exposed within a 600 m thick section of arc-related mantle harzburgites and dunites in the Herbeira massif of the Cabo Ortegal Complex, Spain. We report a combination of field and petrographic observations with in situ and whole-rock geochemical studies of these pyroxenites. After constraining the effects of secondary processes (serpentinization, fluid or melt percolation and amphibolitization), we determine that the low Al content of pyroxenes, high abundance of compatible elements and the absence of plagioclase reflect melt–peridotite interaction and crystal segregation from primitive hydrous melts at relatively low pressure (〈1·2 GPa). Olivine clinopyroxenites and olivine websterites preserving dunite lenses (type 1 and 3 pyroxenites) represent the products of partial replacement of peridotites at decreasing melt/rock ratio following the intrusion of picritic melts. Massive websterites (type 2) may represent the final products of this reaction at higher melt/rock ratios. They crystallized from more Si-rich (boninitic) melts, potentially generated through differentiation of the initially picritic melts or intruded as dykes and veins. Rare opx-rich websterites (type 4) were produced by interaction of these melts with dunites. Chromatographic re-equilibration accompanied late-magmatic crystallization of amphibole from migrating or trapped residual melts. This percolative fractional crystallization produced a range of rare earth element (REE) patterns from spoon-shaped in type 1 pyroxenites to strongly light REE (LREE)-enriched in type 2 and 3 pyroxenites. Particularly high CaO/Al 2 O 3 ratios (2·2–11·3) and the selective enrichment of large ion lithophile elements (LILE) over high field strength elements (HFSE) in Cabo Ortegal pyroxenites suggest the generation of Ca-rich picritic–boninitic parental melts via low-degree, second-stage melting of a refractory lherzolite at 〈2 GPa, following percolation of slab-derived fluids and/or carbonatite melts. Pyroxenites and their host peridotites record high-temperature deformation followed by the development of sheath folds and mylonites. Peak metamorphism was then reached under eclogite-facies conditions (1·6–1·8 GPa and 780–800°C) as recorded by undeformed garnet coronas around spinel. We suggest that this episode corresponds to the delamination of an arc root owing to gravitational instabilities arising from the presence of abundant pyroxenites within mantle harzburgites. Retrograde metamorphism and hydration under amphibolite-facies conditions were recorded by abundant post-kinematic amphibole, which corresponds to the exhumation of the arc root after its intrusion into a subduction zone. The Cabo Ortegal Complex thus preserves a unique section of delaminated arc root, providing evidence for the significant role of melt–peridotite interaction during the differentiation of primitive arc magmas at depth.

Description: Island arc picrites and boninites are magnesian magmatic rocks believed to be generated by high degrees of melting of depleted mantle sources fluxed by subduction-derived, volatile-rich components. These magmas can be probes of both the mantle wedge protoliths and subduction components, but are rare among other, usually more evolved, types of arc lavas. Furthermore, many arc picrites and boninites show evidence for late-stage differentiation prior to or during eruption, masking their primary, mantle-derived geochemical signatures. We report textural and chemical data on spinel-hosted melt inclusions of mantle origin in amphibole-bearing websterite veins cross-cutting spinel harzburgite xenoliths from the active andesitic Avacha volcano (south Kamchatka, Russia). The data are used to constrain the composition and origin of melts that formed the websterite veins in the sub-arc lithospheric mantle. The melt inclusions typically contain euhedral orthopyroxene and clinopyroxene and occasionally minor amphibole in silicate glass. The melt inclusions were homogenized using heating stages and gas-mixing furnaces. The homogenized glasses range from subalkaline primitive silica-rich picrite and high-Ca boninite (〉15 wt % MgO, 48–54 wt % SiO 2 ) to rhyolite. High-Ca boninite glasses have moderate volatile and low heavy rare earth element contents and elevated Cs, Rb, Ba, U, Sr, and Li abundances, with extremely high U/Th. In turn, the glasses display no negative spikes in the high field strength elements Nb, Ta, Zr, Hf, and Ti. We show that the silica-rich picrite and high-Ca boninite liquids in this study formed by high degrees of melting (〉25%), at volatile under-saturation, of hybrid melt-depleted but silica-rich mantle sources at ≥1·5 GPa. The hybrid sources formed in two stages: first, by extraction of ~15% melt from the convecting mantle to form a refractory protolith, which was subsequently enriched in silica via interaction with subduction-derived components prior to or during remelting in the mantle wedge. The subduction-derived components were enriched in fluid-mobile elements and probably oxidized. Overall, our results suggest that silica-rich picrites and high-Ca boninites can be primary melts in mature subduction zones and differentiate within the mantle wedge and the deep arc crust to form more evolved andesite magmas.

Description: Andesitic arc volcanism is the most common type of subduction-related magmatism on Earth. How these melts are generated and under which conditions they evolve towards silica-rich liquids is still a matter of discussion. We have performed crystallization experiments on a representative andesite sample from the Upper Scoriae 1 (USC-1) eruption of Santorini (Greece) with the aim of understanding such processes. Experiments were performed between 1000 and 900 °C, in the pressure range 100–400 MPa, at f O 2 from QFM (quartz–fayalite–magnetite) to NNO (nickel–nickel oxide) + 1·5, with H 2 O melt contents varying from saturation to nominally dry conditions. The results show that the USC-1 andesitic magma was generated at 1000 °C and 12–15 km depth (400 MPa), migrated to shallower levels (8 km; 200 MPa) and intruded into a partially crystallized dacitic magma body. The magma cooled to 975 °C and generated the phenocryst assemblage and compositional zonations that characterize the products of this eruption. An injection of basaltic magma problably subsequently triggered the eruption. In addition to providing the pre-eruptive conditions of the USC-1 magma, our experiments also shed light on the generation conditions of silica-rich magmas at Santorini. Experimental of runs performed at f O 2 ~ NNO + 1 (± 0·5) closely mimic the compositional evolution of magmas at Santorini whereas those at reduced conditions (QFM) do not. Glasses from runs at 1000–975 °C encompass the magma compositions of intermediate-dominated eruptions, whereas those at 950–900 °C reproduce the silicic-dominated eruptions. Altogether, the comparison between our experimental results and natural data for major recent eruptions from Santorini shows that different magma reservoirs, located at different levels, were involved during highly energetic events. Our results suggest that fractionation in deep reservoirs may give rise to magma series with a tholeiitic signature whereas at shallow levels calc-alkaline trends are produced.

Description: Understanding the processes responsible for complex igneous flow structures and layering along the contacts between homogeneous plutons is fundamental for constraining pluton growth and melt transport mechanisms. Orbicules and comb layers from Fisher Lake, Northern Sierra Nevada, are associated with Cretaceous gabbroic to dioritic intrusions. Orbicule cores range from hornblende-bearing olivine-norites, hornblende-norites and hornblende-gabbros/diorites to broken-off comb layers and occasional metamorphic rocks. Magmatic orbicule cores represent cumulates or crystal mush formed by the fractionation of hydrous basaltic andesite magma at upper crustal pressures. We provide field evidence of dykes containing orbicules and/or comb layers as well as magmatic breccias coexisting with orbicules. Orbicule cores, orbicule matrices and host-rocks have overlapping mineral compositions and ages (110–113 Ma), providing evidence for the formation of comb layers and orbicule bodies directly from rising, hydrous basaltic andesite magma. The hornblende-gabbro cumulate signature of the orbicule matrix indicates that the residual melt expelled after the formation of these orbicule bodies had an andesitic composition. We propose, on the basis of field, petrological and geochemical constraints, that the Fisher Lake orbicule bodies were formed through the continuous tapping of a near-homogeneous magma reservoir remobilizing pre-existing, low-pressure, hydrous cumulates prior to injection as orbicule- and diatreme-like breccia dykes at the margin of cooling upper crustal plutons. The Fisher Lake orbicules and comb layers represent subvertical melt extraction zones forming in the roots of a volcanic plumbing system; they are potentially magmatic features linking deep-seated magmatic systems with shallow volcanic environments.

Description: In finite-difference (FD) method, numerical dispersion is the dominant factor influencing the accuracy of seismic modelling. Various optimized FD schemes for scalar wave modelling have been proposed to reduce grid dispersion, while the optimized time–space domain FD schemes for elastic wave modelling have not been fully investigated yet. In this paper, an optimized FD scheme with Equivalent Staggered Grid (ESG) for elastic modelling has been developed. We start from the constant P - and S -wave speed elastic wave equations and then deduce analytical plane wave solutions in the wavenumber domain with eigenvalue decomposition method. Based on the elastic plane wave solutions, three new time–space domain dispersion relations of ESG elastic modelling are obtained, which are represented by three equations corresponding to P -, S - and converted-wave terms in the elastic equations, respectively. By using these new relations, we can study the dispersion errors of different spatial FD terms independently. The dispersion analysis showed that different spatial FD terms have different errors. It is therefore suggested that different FD coefficients to be used to approximate the three spatial derivative terms. In addition, the relative dispersion error in L2-norm is minimized through optimizing FD coefficients using Newton's method. Synthetic examples have demonstrated that this new optimal FD schemes have superior accuracy for elastic wave modelling compared to Taylor-series expansion and optimized space domain FD schemes.

Description: Markov chain Monte-Carlo (McMC) sampling generates correlated random samples such that their distribution would converge to the true distribution only as the number of samples tends to infinity. In practice, McMC is found to be slow to converge, convergence is not guaranteed to be achieved in finite time, and detection of convergence requires the use of subjective criteria. Although McMC has been used for decades as the algorithm of choice for inference in complex probability distributions, there is a need to seek alternative approaches, particularly in high dimensional problems. Walker & Curtis ( 2014 ) developed a method for Bayesian inversion of 2-D spatial data using an exact sampling alternative to McMC which always draws independent samples of the target distribution. Their method thus obviates the need for convergence and removes the concomitant bias exhibited by finite sample sets. Their algorithm is nevertheless computationally intensive and requires large memory. We propose a more efficient method for Bayesian inversion of categorical variables, such as geological facies that requires no sampling at all. The method is based on a 2-D Hidden Markov Model (2D-HMM) over a grid of cells where observations represent localized data constraining each cell. The data in our example application are seismic attributes such as P - and S -wave impedances and rock density; our categorical variables are the hidden states and represent the geological rock types in each cell—facies of distinct subsets of lithology and fluid combinations such as shale, brine-sand and gas-sand. The observations at each location are assumed to be generated from a random function of the hidden state (facies) at that location, and to be distributed according to a certain probability distribution that is independent of hidden states at other locations – an assumption referred to as ‘localized likelihoods’. The hidden state (facies) at a location cannot be determined solely by the observation at that location as it also depends on prior information concerning the spatial distribution of other hidden states elsewhere. The prior information is included in the inversion in the form of a training image which represents a conceptual depiction of local geologies that might be expected, but other forms of prior information can be used in the method as desired. The method provides direct estimates of posterior marginal probability distributions over each variable, so these do not need to be estimated from samples such as in McMC. Nevertheless, in the case that samples are desired, these can be generated. On a 2-D test example the method is shown to outperform previous methods significantly, at a fraction of the computational cost. In many foreseeable applications there are therefore no serious impediments to extending the method to 3-D cases.

Description: Using an up-to-date global plate rotation model, applied to the endpoints of preserved major spreading ridge isochrons, we have calculated the explicitly reconstructable length-weighted mean global half-spreading rate (HSR), ridge length and area production as a function of time since the end of the Cretaceous Normal Superchron at 83.0 Ma. Our calculations integrate uncertainties in rotation parameters and chron boundary ages with the partial sampling uncertainties arising from progressive subduction of older oceanic lithosphere and its preserved spreading record. This record of directly reconstructable oceanic ridge production provides a well-constrained baseline that can be compared to reconstructions that include the largely unconstrained extrapolated histories of entirely subducted oceanic plates. The directly reconstructable global mean HSR has not varied by more than ±15 per cent about an average rate of 28.4 ± 4.6 mm a –1 since 83 Ma. No long-term secular trend is evident: a maximum global mean half-rate of 32 ± 6 mm a –1 occurred from 33.1 Ma to about 25.8 Ma, with minima of 26 ± 5 mm a –1 between about 56 and 40.2 Ma, and 24 ± 1 mm a –1 since 3.2 Ma. Only this most recent interval has a rate that differs significantly (at ±2) from the long-term mean. The global, reconstructable ridge length at 56 Ma decreases by less than 15 per cent relative to the modern ridge system; by 83 Ma it has decreased by 38 per cent. These relatively high preserved ridge fractions mean that the estimated uncertainty due to partial sampling stays roughly equivalent to the estimated rotation model uncertainties, allowing long-term spreading rate variations of 〉20 per cent since the Late Cretaceous to be ruled out. In contrast, prior to 83 Ma too little oceanic lithosphere is preserved to reliably reconstruct global spreading rates.

Description: A double-correlation method is introduced to locate tremor sources based on stacks of complex, doubly-correlated tremor records of multiple triplets of seismographs back projected to hypothetical source locations in a geographic grid. Peaks in the resulting stack of moduli are inferred source locations. The stack of the moduli is a robust measure of energy radiated from a point source or point sources even when the velocity information is imprecise. Application to real data shows how double correlation focuses the source mapping compared to the common single correlation approach. Synthetic tests demonstrate the robustness of the method and its resolution limitations which are controlled by the station geometry, the finite frequency of the signal, the quality of the used velocity information and noise level. Both random noise and signal or noise correlated at time shifts that are inconsistent with the assumed velocity structure can be effectively suppressed. Assuming a surface wave velocity, we can constrain the source location even if the surface wave component does not dominate. The method can also in principle be used with body waves in 3-D, although this requires more data and seismographs placed near the source for depth resolution.

Description: The possibility of applying one explicit finite-difference (FD) scheme to all interior grid points (points not lying on a grid border) no matter what their positions are with respect to the material interface is one of the key factors of the computational efficiency of the FD modelling. Smooth or discontinuous heterogeneity of the medium is accounted for only by values of the effective grid moduli and densities. Accuracy of modelling thus very much depends on how these effective grid parameters are evaluated. We present an orthorhombic representation of a heterogeneous medium for the FD modelling. We numerically demonstrate its superior accuracy. Compared to the harmonic-averaging representation the orthorhombic representation is more accurate mainly in the case of strong surface waves that are especially important in local surface sedimentary basins. The orthorhombic representation is applicable to modelling seismic wave propagation and earthquake motion in isotropic models with material interfaces and smooth heterogeneities using velocity–stress, displacement–stress and displacement FD schemes on staggered, partly staggered, Lebedev and collocated grids.

Description: We study fluctuations in the degree-2 zonal spherical harmonic coefficient of the Earth's gravity potential, C 20 , over the period 2003–2015. This coefficient is related to the Earth's oblateness and studying its temporal variations, C 20 , can be used to monitor large-scale mass movements between high and low latitude regions. We examine C 20 inferred from six different sources, including satellite laser ranging (SLR), GRACE and global geophysical fluids models. We further include estimates that we derive from measured variations in the length-of-day (LOD), from the inversion of global crustal displacements as measured by GPS, as well as from the combination of GRACE and the output of an ocean model as described by Sun et al. We apply a sequence of trend and seasonal moving average filters to the different time-series in order to decompose them into an interannual, a seasonal and an intraseasonal component. We then perform a comparison analysis for each component, and we further estimate the noise level contained in the different series using an extended version of the three-cornered-hat method. For the seasonal component, we generally obtain a very good agreement between the different sources, and except for the LOD-derived series, we find that over 90 per cent of the variance in the seasonal components can be explained by the sum of an annual and semiannual oscillation of constant amplitudes and phases, indicating that the seasonal pattern is stable over the considered time period. High consistency between the different estimates is also observed for the intraseasonal component, except for the solution from GRACE, which is known to be affected by a strong tide-like alias with a period of about 161 d. Estimated interannual components from the different sources are generally in agreement with each other, although estimates from GRACE and LOD present some discrepancies. Slight deviations are further observed for the estimate from the geophysical models, likely to be related to the omission of polar ice and groundwater changes in the model combination we use. On the other hand, these processes do not seem to play an important role at seasonal and shorter timescales, as the sum of modelled atmospheric, oceanic and hydrological effects effectively explains the observed C 20 variations at those scales. We generally obtain very good results for the solution from SLR, and we confirm that this well-established technique accurately tracks changes in C 20 . Good agreement is further observed for the estimate from the GPS inversion, showing that this indirect method is successful in capturing fluctuations in C 20 on scales ranging from intra- to interannual. Obtaining accurate estimates from LOD, however, remains a challenging task and more reliable models of atmospheric wind fields are needed in order to obtain high-quality C 20 , in particular at the seasonal scale. The combination of GRACE data and the output of an ocean model appears to be a promising approach, particularly since corresponding C 20 is not affected by tide-like aliases, and generally gives better results than the solution from GRACE, which still seems to be of rather poor quality.

Description: The diffusion of titanium along with Ni, Co and Mn in pure synthetic forsterite has been studied as a function of temperature, oxygen fugacity, crystallographic orientation and chemical potentials in the four-component system MgO–SiO 2 –TiO 2 –TiO 1·5 . In over 100 different experimental conditions, no Ti diffusion profiles with the classical shape of the error function were observed. Instead, the profiles vary between ‘stepped’ at high f O 2 (probably owing to diffusion on octahedral sites and trapping in tetrahedral sites), to hockey-stick shapes at low f O 2 owing to concentration-dependent diffusion. The profile shapes vary systematically between these end-members according to oxygen fugacity. These profile shapes are also reproduced in experiments using natural San Carlos olivine. The change in profile shape is interpreted to be a function of valence-state change from Ti 4+ to Ti 3+ , shown by a sigmoidal log f O 2 –diffusivity relationship. To distinguish between Ti 4+ and Ti 3+ in forsterite, a diffusion profile was also investigated using ‘hydroxylation spectroscopy’ whereby a crystal from a diffusion experiment was hydroxylated by annealing in an H 2 O-rich fluid at moderate pressure and temperature. The point defects associated with the added structural OH were then determined by infrared spectroscopy, providing evidence for the presence of Ti 4+ and Ti 3+ at intermediate oxygen fugacity. The transition from Ti 4+ to Ti 3+ occurs at considerably higher f O 2 in these diffusion experiments than in equilibrium experiments in similar systems. Titanium diffusion, qualified using a mobility parameter ( M ), is fastest along the c -axis, at high activity of silica, low oxygen fugacity and high temperature. The rate of Ti diffusion is broadly similar to those of Mn, Ni and Co, and closer to published rates of Mg self-diffusion in olivine than Si self-diffusion. The observed Ti 3+ and Ti 4+ diffusion occurs on the M-sites. From statistical examination of the large Mn, Ni and Co diffusivity dataset, we determine that there is a 10 2/3 dependence of diffusion on a SiO 2 ; the cations diffuse more rapidly at high a SiO 2 . In addition, diffusivity has a 10 –1/8 dependence on f O 2 , suggested to be the result of increased concentration of loosely bound Ti 3+ –vacancy pairs, which enhance the mobility of other cations. This behaviour is likely to be present in natural systems where trivalent (e.g. Fe 3+ , Al 3+ , Cr 3+ ) and divalent cations diffuse together. Consideration should be given to diffusive interference by highly charged, fast-moving cations when extracting time scales from frozen diffusion profiles.

Description: The role of lithospheric mantle metasomatized by CO 2 -bearing melts in the genesis of HIMU-like alkaline intraplate basalts is investigated using a suite of peridotite xenoliths from New Zealand. The xenoliths have Sr–Nd–Pb–Hf isotope compositions ( 87 Sr/ 86 Sr = 0·7029, Nd = + 5 to + 6, 206 Pb/ 204 Pb = 20·4 and Hf = +5 to +8) indistinguishable from the host low-silica basalts and, except for 207 Pb/ 204 Pb, overlapping with the HIMU mantle reservoir. Laser line scans across grain boundaries in the xenoliths show, however, that the host magma contribution is restricted to minor degrees of melt infiltration along grain boundaries during ascent, with the distinctive peridotite isotopic compositions having been imparted earlier by mantle metasomatism. Two mantle metasomatic styles are distinguished from pyroxene trace element concentrations (in particular, rare earth elements, Ti, Zr and Hf) and are interpreted to be the result of reaction of peridotite with CO 2 - bearing magmas. The occurrence of two subtly chemically different but isotopically indistinguishable styles of metasomatism in rocks with the same equilibrium temperatures within the same mantle column may be due to separate volatile-rich melts formed by different degrees of melting of a deeper carbonated peridotitic ± pyroxenitic source, or due to metasomatism having been imparted to different degrees on a variably depleted protolith. In either case, the formation of the HIMU-like enriched lithospheric mantle was achieved by percolation of volatile-rich melts, which probably rose from the asthenosphere. Melt modelling of representative depleted and subsequently enriched samples shows that low-degree melting of a CO 2 -bearing melt-metasomatized peridotite could yield a melt with a trace element composition very similar to that of the Zealandia HIMU-like alkaline basalts, but only if small volumes (~5%) of amphibole participated in the melting process. Although not observed in the studied xenoliths, amphibole is associated with mantle metasomatism by carbonatitic or CO 2 -bearing melts elsewhere in the world and has been found as xenocrysts with HIMU-like isotope compositions in some Zealandia basalts. The melt modelling results also imply that amphibole could buffer the trace element budgets of a low-degree melt regardless of the source peridotite composition; therefore, provided that hydrous metasomatized lithospheric mantle can be perturbed to melt, the contribution of amphibole would explain the similarities of alkaline ocean island basalt-like magmas in continental and oceanic settings. HIMU-like reservoirs formed by percolation of young volatile-rich (CO 2 + H 2 O) melts are widespread within the Earth’s lithospheric mantle and are a potential source for intraplate alkaline basaltic magmatism.

Description: Joint inversion that simultaneously inverts multiple geophysical data sets to recover a common Earth model is increasingly being applied to exploration problems. Petrophysical data can serve as an effective constraint to link different physical property models in such inversions. There are two challenges, among others, associated with the petrophysical approach to joint inversion. One is related to the multimodality of petrophysical data because there often exist more than one relationship between different physical properties in a region of study. The other challenge arises from the fact that petrophysical relationships have different characteristics and can exhibit point, linear, quadratic, or exponential forms in a crossplot. The fuzzy c-means (FCM) clustering technique is effective in tackling the first challenge and has been applied successfully. We focus on the second challenge in this paper and develop a joint inversion method based on variations of the FCM clustering technique. To account for the specific shapes of petrophysical relationships, we introduce several different fuzzy clustering algorithms that are capable of handling different shapes of petrophysical relationships. We present two synthetic and one field data examples and demonstrate that, by choosing appropriate distance measures for the clustering component in the joint inversion algorithm, the proposed joint inversion method provides an effective means of handling common petrophysical situations we encounter in practice. The jointly inverted models have both enhanced structural similarity and increased petrophysical correlation, and better represent the subsurface in the spatial domain and the parameter domain of physical properties.

Description: We present two independent, automated methods for estimating the absolute horizontal misorientation of seismic sensors. We apply both methods to 44 free-fall ocean-bottom seismometers (OBSs) of the RHUM-RUM experiment ( http://www.rhum-rum.net/ ). The techniques measure the 3-D directions of particle motion of (1) P -waves and (2) Rayleigh waves of earthquake recordings. For P -waves, we used a principal component analysis to determine the directions of particle motions (polarizations) in multiple frequency passbands. We correct for polarization deviations due to seismic anisotropy and dipping discontinuities using a simple fit equation, which yields significantly more accurate OBS orientations. For Rayleigh waves, we evaluated the degree of elliptical polarization in the vertical plane in the time and frequency domain. The results obtained for the RHUM-RUM OBS stations differed, on average, by 3.1° and 3.7° between the methods, using circular mean and median statistics, which is within the methods’ estimate uncertainties. Using P -waves, we obtained orientation estimates for 31 ocean-bottom seismometers with an average uncertainty (95 per cent confidence interval) of 11° per station. For 7 of these OBS, data coverage was sufficient to correct polarization measurements for underlying seismic anisotropy and dipping discontinuities, improving their average orientation uncertainty from 11° to 6° per station. Using Rayleigh waves, we obtained misorientation estimates for 40 OBS, with an average uncertainty of 16° per station. The good agreement of results obtained using the two methods indicates that they should also be useful for detecting misorientations of terrestrial seismic stations.

Description: It is shown that a slip wave solution exists for antiplane sliding of an elastic layer on an elastic half-space. It is a companion solution to the well-known Love wave solution.

Description: Surface waves propagating in layered media inherently possess multimodal dispersion characteristics. However, traditional surface wave testing methods employing measurements at the free surface usually capture only a single apparent dispersion curve, especially when using short geophone arrays common to near surface and geotechnical-scale investigations. Such single-mode or fragmented multimode apparent dispersion curves contain only a fraction of the possible dispersion information, thus limiting the accuracy of inverted profiles. To enable more robust measurement of higher Rayleigh-wave modes, a recently developed hybrid minimally invasive multimodal surface wave method is combined herein with the widely used geotechnical standard penetration test (SPT), which is employed as a practical and ubiquitous downhole source. Upon superimposing surface wave dispersion data for a range of SPT impact depths within the soil, higher modes can be measured more consistently and reliably relative to traditional non-invasive testing methods. As a result, misidentification of multiple dispersion modes can be practically eliminated, significantly improving the accuracy and certainty of inversion results.

Description: Coal seam fires are a worldwide disaster of both ecological and economic importance. Their remote detection from the ground surface or using airborne techniques is required for developing efficient strategies to extinguish them. We investigate here the use of time-domain-induced polarization to localize coal seam fires. For laboratory experiments, we first introduce a modified time-domain-induced polarization methodology to quickly acquire and invert the secondary voltage distribution mapped after the shutdown of the primary current. A set of sandbox experiments is conducted in which coal is embedded into humidified sand. Raw coal alone generates significant induced polarization anomalies, above those shown by the sand. Even higher induced polarization anomalies are detected in presence of a coal seam fire. We postulate that the higher chargeability is due to the pyrolysis, which may enhance electronic polarization or the polarization associated with the cation exchange capacity (CEC) of the material. The position of the coal seam fire is well recovered inside the tank by inverting the secondary voltages in term of a source current density distribution. We also collected field data over a recognized coal seam fire in Colorado, USA. A chargeability anomaly (~800 mV V –1 ) and a resistivity anomaly (~1 Ohm m) are observed at the position of the coal seam fire. We propose a normalized burning front index (a scaled normalized chargeability) to image and localize, without ambiguity, the position of the coal seam fire in the subsurface. The 3-D reconstructed target is located below a negative self-potential anomaly (similarly to what is observed in laboratory experiments) and a temperature anomaly recorded at a depth of 30 cm.

Description: We investigate the elastic and anelastic structure of the lowermost mantle at the western edge of the Pacific large low shear velocity province (LLSVP) by inverting a collection of S and ScS waveforms. The transverse component data were obtained from F-net for 31 deep earthquakes beneath Tonga and Fiji, filtered between 12.5 and 200 s. We observe a regional variation of S and ScS arrival times and amplitude ratios, according to which we divide our region of interest into three subregions. For each of these subregions, we then perform 1-D (depth-dependent) waveform inversions simultaneously for radial profiles of shear wave velocity ( V S ) and seismic quality factor ( Q ). Models for all three subregions show low V S and low Q structures from 2000 km depth down to the core–mantle boundary. We further find that V S and Q in the central subregion, sampling the Caroline plume, are substantially lower than in the surrounding regions, whatever the depth. In the central subregion, V S -anomalies with respect to PREM (d V S ) and Q are about –2.5 per cent and 216 at a depth of 2850 km, and –0.6 per cent and 263 at a depth of 2000 km. By contrast, in the two other regions, d V S and Q are –2.2 per cent and 261 at a depth of 2850 km, and –0.3 per cent and 291 at a depth of 2000 km. At depths greater than ~2500 km, these differences may indicate lateral variations in temperature of ~100 K within the Pacific LLSVP. At shallower depths, they may be due to the temperature difference between the Caroline plume and its surroundings, and possibly to a small fraction of iron-rich material entrained by the plume.