ALBERT

Description: 〈span〉〈div〉Summary〈/div〉The strength of the lithosphere plays a key role in the formation and evolution of tectonic plate boundaries. Localized lithospheric deformation associated with plate tectonics requires a mechanism for weakening across the entire width of the lithosphere, including the strongest cold ductile region. We explore the microphysics of weakening of lithospheric materials, and in particular the coupled evolution of mineral grain size and intragranular defects and their control on lithospheric strength. We propose a model for the interaction between grain-boundaries and dislocation density to reduce the net free energy of grains during dynamic recrystallization (DRX). The driving forces for DRX arise from heterogeneity in dislocation density and grain boundary curvature. Our model shows that grain growth driven by variation in grain boundary curvature can be impeded by variation in dislocation density; this occurs because as the grains grow, to minimize their surface energy, their dislocation density and associated internal energy may increase and offset the driving forces for grain growth. The correlation between grain size and dislocation density can for example arise because the dislocation accumulation in smaller grains is suppressed due to the large stress that is needed to bend and elongate a short dislocation (as dictated by the small grain size), while the larger grains can have long dislocations and reach a steady state dislocation density dictated by the applied stress. In a lithospheric setting, slower grain growth means that it would require less mechanical work to establish weak localized shear zones through grain damage, and retard the healing of previously damaged zones. Furthermore, the competition of two different time-scales - that of grain growth and the dislocation kinetics - can lead to oscillating behavior over 1 to 10 years as the grain size and dislocation density advance towards their steady states. These oscillations are likely to have an effect on the rheology of lithospheric rocks, e.g. their strengthening and weakening through time, and have a potential application to geological processes such as postseismic creep in ductile shear zones.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Locating and monitoring passive seismic sources provides us important information for studying subsurface rock deformation, injected fluid migration, regional stress conditions as well as fault rupture mechanism. In this paper, we present a novel passive-source monitoring approach using vector-based elastic time-reversal imaging. By solving the elastic wave equation using observed multicomponent records as boundary conditions, we first compute back-propagated elastic wavefields in the subsurface. Then, we separate the extrapolated wavefields into compressional (P-wave) and shear (S-wave) modes using the vector Helmholtz decomposition. A zero-lag cross-correlation imaging condition is applied to the separated pure-mode vector wavefields to produce passive-source images. We compare imaging results using three implementations, i.e., dot-product, energy and power. Numerical experiments demonstrate that the power imaging condition gives us the highest resolution and is less sensitive to the presence of random noises. To capture the propagation of microseismic fracture and earthquake rupture, we modify the traditional zero-lag cross-correlation imaging condition by summing the multiplication of the separated P- and S-wavefields within local time windows, which enables us to capture the temporal and spatial evolution of earthquake rupture. 2D and 3D numerical examples demonstrate that the proposed method is capable of accurately locating point sources, as well as delineating dynamic propagation of hydraulic fracture and earthquake rupture.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Magnitudes of differential stress in the lithosphere, especially in the crust, are still disputed. Earthquake-based stress drop estimates indicate median values 〈 10 MPa whereas the lateral variation of gravitational potential energy per unit area (〈span〉GPE〈/span〉) across significant relief indicates stress magnitudes of ca. 100 MPa in average across a 100 km thick lithosphere between the Indian lowland and the Tibetan plateau. These standard 〈span〉GPE〈/span〉-based stress estimates correspond to membrane stresses, because they are associated with a deformation that is uniform with depth. We show here with new analytical results that lateral variations in 〈span〉GPE〈/span〉 can also cause bending moments and related bending stresses of several hundreds of MPa. Furthermore, we perform two-dimensional thermo-mechanical numerical simulations (1) to evaluate estimates for membrane and bending stresses based on 〈span〉GPE〈/span〉 variations, (2) to quantify minimum crustal stress magnitudes that are required to maintain the topographic relief between Indian lowland and Tibetan plateau for ca. 10 Ma and (3) to quantify the corresponding relative contribution of crustal strength to the total lithospheric strength. The numerical model includes viscoelastoplastic deformation, gravity and heat transfer. The model configuration is based on density fields from the CRUST1.0 data set and from a geophysically and petrologically constrained density model based on 〈span〉in situ〈/span〉 field campaigns. The numerical results indicate that values of differential stress in the upper crust must be 〉 ca. 180 MPa, corresponding to a friction angle of ca. 10°, to maintain the topographic relief between lowland and plateau for 〉 10 Ma. The relative contribution of crustal strength to total lithospheric strength varies considerably laterally. In the region between lowland and plateau and inside the plateau the depth-integrated crustal strength is approximately equal to the depth-integrated strength of the mantle lithosphere. Simple analytical formulae predicting the lateral variation of depth-integrated stresses agree with numerically calculated stress fields, which show both the accuracy of the numerical results and the applicability of simple, rheology-independent, analytical predictions to highly variable, rheology-dependent, stress fields. Our results indicate that (1) crustal strength can be locally equal to mantle lithosphere strength and that (2) crustal stresses must be at least one order of magnitude larger than median stress drops in order to support the plateau relief over a duration of ca. 10 Ma.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Scanning magnetometers are increasingly used to characterize the magnetization of mineral grains in rock samples. Up-scaling this measurement technique to large numbers of individual particles is hampered by the intrinsic non-uniqueness of potential-field inversion. Here it is shown that this problem can be circumvented by adding tomographic information that determines the location of the possible field sources. Standard potential theory is used to prove a uniqueness theorem which completely characterizes the mathematical background of the corresponding source-localized inversion. It exactly resolves under which conditions a potential field measurement on a surface can be uniquely decomposed into signals from the different source regions. The intrinsic non-uniqueness of potential field inversion prevents that the source distribution inside the tomographically outlined regions can be recovered, but the potential field of each region is uniquely defined. For scanning magnetometers in rock magnetism, this result implies that magnetic dipole vectors of large numbers of individual magnetic particles can be reliably reconstructed from surface scans of the magnetic field, if the particle positions are independently determined. This provides an incentive to improve scanning methods for future paleomagnetic applications.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The most common earthquake forecasting models assume that the magnitude of the next earthquake is independent from the past. This feature severely limits the capability to forecast large earthquakes with high probabilities. Here we investigate empirically on the magnitude-independence assumption, exploring if: i) background and triggered earthquakes have the same frequency-magnitude distribution, ii) variations of seismicity in the space-time-magnitude domain encode some information on the future earthquakes size. For this purpose, and to verify the stability of the findings, we consider seismic catalogues covering different space-time-magnitude windows, such as the Alto Tiberina Near Fault Observatory (TABOO), the California and Japanese seismic catalogues. Our approach is inspired by the nearest-neighbour method proposed by Baiesi & Paczuski (2004) and elaborated by Zaliapin et al. (2008) to distinguish between triggered and background earthquakes. Here we implement the same metric-based correlation to identify the precursory seismicity of any triggered earthquake; this allows us to analyse, for each triggered earthquake, the space-time-magnitude distribution of the seismicity that likely contributed to its occurrence. Our results show that the magnitude-independence assumption holds reasonably well in all catalogues, with a remarkable exception that is consistent with a previous independent study; this departure from the magnitude-independence assumption shows that larger events tend to nucleate at a higher distance from the ongoing sequence. We also notice that the reliability of this assumption may depend on the spatial scale considered; it holds for seismic catalogues of large areas, but we identify possible departures in small areas, reflecting different ways to release locally seismic energy. Finally, we come across an important issue that may lead to misleading results in similar studies, i.e., if a seismic catalogue appears overall complete above a fixed magnitude threshold, it may still yield spurious signals into the analysis. Specifically, we show that some significant departures from the magnitude-independence assumption do not survive when considering spatiotemporal variations of the magnitude of completeness.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Vp/Vs models provide important complementary information to Vp and Vs models, relevant to lithology, rock damage, partial melting, water saturation, etc. However, seismic tomography using body-wave traveltime data from local or regional earthquakes does not constrain Vp/Vs well due to the different resolution of Vp and Vs models, with the Vp models usually better constrained than Vs. Since surface-wave data are most sensitive to Vs, which leads to complementary strengths with respect to body-wave data, we jointly invert body- and surface-wave data to better resolve the Vp/Vs models. In order to show the robustness of our joint inversion method, we compare the results to other approaches, including dividing Vp by Vs models and Vp/Vs parameterization with body-wave or both body- and surface-wave data, using synthetic data and real data from the southern California plate boundary region. We confirm that Vp/Vs models from body-wave inversion obtained by division tend to include artifacts, even when both Vp and Vs models seem reasonable. The joint inversion with Vp/Vs parameterization is found to improve the Vp/Vs ratio model significantly and the resultant Vp/Vs model shows more geologically consistent features, such as high Vp/Vs along fault traces at shallow depths likely indicating fault-related damage. The Vp/Vs model also exhibits contrasts at intermediate depths along the Peninsular Range compositional boundary, and high Vp/Vs in the lower crust near the Salton Sea region correlated with high heat flow and may indicate partial melting. The improved Vp/Vs as well as individual Vp and Vs models are useful for earthquake relocation, high-resolution Moho depth imaging, and interpretation of other data and tectonic evolution in the region.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Interpretation of surface fault scarps and palaeoseismic trenches is a key component of estimating fault slip rates, earthquake recurrence rates and maximum magnitudes for hazard assessments. Often these analyses rely on the assumption that successive earthquakes all breached the surface and that the ruptures are recorded topographically, or by the deposits exposed in a trench. The M〈sub〉〈span〉w〈/span〉〈/sub〉7.2 1992 Suusamyr earthquake, Kyrgyzstan, is an apparently problematic case for such analyses because its ruptures show significant displacement but are only mapped as having broken the surface along small, disparate portions of the fault. Here we present the results of surveys conducted along the Suusamyr Fault to establish whether that is the case. Two sets of ruptures were identified following the earthquake. They are unusually short for their displacement and are separated by a 25 km gap. Using satellite imagery, high-resolution digital elevation models and palaeoseismic trenching we first reassess the distribution of the 1992 ruptures and then reconstruct the Holocene earthquake record to establish the extent to which the 1992 earthquake is representative of the rupture behaviour of this fault. We find evidence for at least two prehistoric surface rupturing earthquakes in the Holocene: one ∼3 ka and one 〉8 ka that, along with the modern event, gives recurrence intervals of ∼3 kyr and ∼5 kyr. Within spatial gaps between segments of the 1992 ruptures there are clear prehistoric surface ruptures and the ruptures in each prehistoric earthquake were discontinuous. We conclude that there is significant variability in the surface rupture pattern of successive earthquakes on the Suusamyr Fault, with implications for the completeness of palaeoseismic records obtained from thrust scarps.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We provide a six-component (6-C) polarization model for 〈span〉P〈/span〉-, 〈span〉SV〈/span〉-, 〈span〉SH〈/span〉-, Rayleigh-, and Love-waves both inside an elastic medium as well as at the free surface. It is shown that single-station 6-C data comprised of three components of rotational motion and three components of translational motion provide the opportunity to unambiguously identify the wave type, propagation direction, and local 〈span〉P〈/span〉- and 〈span〉S〈/span〉-wave velocities at the receiver location by use of polarization analysis. To extract such information by conventional processing of three-component (3-C) translational data would require large and dense receiver arrays. The additional rotational components allow the extension of the rank of the coherency matrix used for polarization analysis. This enables us to accurately determine the wave type and wave parameters (propagation direction and velocity) of seismic phases, even if more than one wave is present in the analysis time window. This is not possible with standard, pure-translational 3-C recordings. In order to identify modes of vibration and to extract the accompanying wave parameters, we adapt the multiple signal classification algorithm (MUSIC). Due to the strong nonlinearity of the MUSIC estimator function, it can be used to detect the presence of specific wave types within the analysis time window at very high resolution. We show how the extracted wavefield properties can be used, in a fully automated way, to separate the wavefield into its different wave modes using only a single 6-C recording station. As an example, we apply the method to remove surface wave energy while preserving the underlying reflection signal and to suppress energy originating from undesired directions, such as side-scattered waves.〈/span〉

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〉The Sentinel-1 mission comprises two synthetic aperture radar satellites, each with a 12 day orbital repeat, orbiting 6 days apart within a narrow tube. The mission design promises the ability to respond quickly to earthquakes with InSAR, and to facilitate production of interferograms with good interferometric correlation globally. We report on our efforts to study global seismicity using Sentinel-1 Interferometric Wide-Swath data between April 2015 and December 2016. We select 35 potentially detectable terrestrial earthquakes in the range 5.5 ≤ 〈span〉Mw〈/span〉 ≤ 7.8 on the basis of their locations, depths and magnitudes, and process the first post-event interferogram with the shortest possible time-span for each using the ISCE software. We evaluate each interferogram for earthquake deformation signals by visual inspection. We can identify deformation signals attributable to earthquakes in 18 of these interferograms (51%); a further six interferograms (17%) have ambiguous interferometric phase affected by tropospheric noise. 11 events (31%) could not be identified from their interferograms. The majority of these failed detections were due to interferogram decorrelation, particularly apparent for earthquakes that occurred between 15°N and 15°S, where climate conditions promote dense vegetation. The majority of the ambiguous interferograms are affected by tropospheric noise, suggesting that techniques to mitigate such noise could improve detection performance. The largest event we do not detect with Sentinel-1 data is a 〈span〉Mw〈/span〉7.0 earthquake that occurred in Vanuatu in April 2016; we also fail to detect the 2016 〈span〉Mw〈/span〉6.2 Kurayoshi earthquake in one out of two possible 24-day interferograms. We propose these as upper and lower estimates on the magnitude of completeness for earthquakes studied with Sentinel-1 data; to raise the magnitude of completeness we suggest that more frequent (e.g. six day) recurrence may be necessary in low latitude areas.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The correct estimation of site-specific attenuation is crucial for the assessment of seismic hazard. Downhole instruments provide in this context valuable information to constrain attenuation directly from data. In this study, we apply an interferometric approach to this problem by deconvolving seismic motions recorded at depth with those recorded at the surface. In doing so, incident and surface-reflected waves can be separated. We apply this technique not only to earthquake data but also to recordings of ambient vibrations. We compute the transfer function between incident and surface-reflected waves in order to infer frequency dependent quality factors for S-waves. The method is applied to a 87 m deep borehole sensor and a co-located surface instrument situated at a hard-rock site in West Bohemia/Vogtland, Germany. We show that the described method provides comparable attenuation estimates using either earthquake data or ambient noise for frequencies between 5-15 Hz. Moreover, a single hour of noise recordings seems to be sufficient to yield stable deconvolution traces and quality factors, thus, offering a fast and easy way to derive attenuation estimates from borehole recordings even in low to mid seismicity regions.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Injection of CO〈sub〉2〈/sub〉 into tight reservoirs produces both gravity change and ground deformation, which provide great opportunities for more accurate coupled inverse modelling. In this study, we incorporate signals generated from several synthetic models to estimate the CO〈sub〉2〈/sub〉 distribution in the reservoir. A relationship is found that connects density variations to volumetric changes associated with injected CO〈sub〉2〈/sub〉, taking advantage of a common set of model parameters for both gravitational and geo-mechanical inverse modelling. This is achieved by assuming that the injected CO〈sub〉2〈/sub〉 increases pressure in the reservoir, which in turn generates extra porosity that is then filled in by the CO〈sub〉2〈/sub〉 mass in the generated space. Tikhonov regularization, supported by the Generalized Cross Validation (GCV) technique for finding the optimized model, is used to solve the ill-posed inverse problems. The results indicate that with a combination of gravity and ground deformation monitoring, the uncertainty and ambiguity in gravimetric modelling due to high levels of noise is mitigated by implementing highly accurate ground deformation measurements, which normally have a higher signal to noise ratio.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉As the number of seismic sensors grows, it is becoming increasingly difficult for analysts to pick seismic phases manually and comprehensively, yet such efforts are fundamental to earthquake monitoring. Despite years of improvements in automatic phase picking, it is difficult to match the performance of experienced analysts. A more subtle issue is that different seismic analysts may pick phases differently, which can introduce bias into earthquake locations. We present a deep-neural-network-based arrival-time picking method called ”PhaseNet” that picks the arrival times of both P and S waves. Deep neural networks have recently made rapid progress in feature learning, and with sufficient training, have achieved super-human performance in many applications. PhaseNet uses three-component seismic waveforms as input and generates probability distributions of P arrivals, S arrivals, and noise as output. We engineer PhaseNet such that peaks in the probability distributions provide accurate arrival times for both P and S waves. PhaseNet is trained on the prodigious available data set provided by analyst-labeled P and S arrival times from the Northern California Earthquake Data Center. The dataset we use contains more than seven hundred thousand waveform samples extracted from over thirty years of earthquake recordings. We demonstrate that PhaseNet achieves much higher picking accuracy and recall rate than existing methods when applied to the waveforms of known earthquakes, which has the potential to increase the number of S-wave observations dramatically over what is currently available. This will enable both improved locations and improved shear wave velocity models.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The low frequency earthquakes (LFEs) that constitute tectonic tremor are often inferred to be slow: to have durations of 0.2 to 0.5 s, a factor of 10 to 100 longer than those of typical 〈span〉MW〈/span〉 1-2 earthquakes. Here we examine LFEs near Parkfield, CA in order to assess several proposed explanations for LFEs’ long durations. We determine LFE rupture areas and location distributions using a new approach, similar to directivity analysis, where we examine how signals coming from various locations within LFEs’ finite rupture extents create differences in the apparent source time functions recorded at various stations. We use synthetic ruptures to determine how much the LFE signals recorded at each station would be modified by spatial variations of the source-station travel time within the rupture area given various possible rupture diameters, and then compare those synthetics with the data. Our synthetics show that the methodology can identify inter-station variations created by heterogeneous slip distributions or complex rupture edges, and thus lets us estimate LFE rupture extents for unilateral or bilateral ruptures. To obtain robust estimates of the sources’ similarity across stations, we stack signals from thousands of LFEs, using an empirical Green’s function approach to isolate the LFEs’ apparent source time functions from the path effects. Our analysis of LFEs in Parkfield implies that LFEs’ apparent source time functions are similar across stations at frequencies up to 8 to 16 Hz, depending on the family. The inter-station coherence observed at these relatively high frequencies, or short wavelengths (down to 0.2 to 0.5 km), suggest that LFEs in each of the 7 families examined occur on asperities. They are clustered in patches with sub-1-km diameters. The individual LFEs’ rupture diameters are estimated to be smaller than 1.1 km for all families, and smaller than 0.5 km and 1 km for the two shallowest families, which were previously found to have 0.2-s durations. Coupling the diameters with the durations suggests that it is possible to model these 〈span〉MW〈/span〉 1-2 LFEs with earthquake-like rupture speeds: around 70% of the shear wave speed. However, that rupture speed matches the data only at the edge of our uncertainty estimates for the family with highest coherence. The data for that family are better matched if LFEs have rupture velocities smaller than 40% of the shear wave speed, or if LFEs have different rupture dynamics. They could have long rise times, contain composite sub-ruptures, or have slip distributions that persist from event to event.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉On June 24, 2015, a 230,000 cubic metre landslide slid into the triangle bayou at the intersection of the Yangtze and Daning Rivers and generated a river tsunami that ran up 6.2 metres on the opposite shoreline at Wushan town. The slope failure and resulting waves killed two people and damaged many shipping facilities. Based on field surveys and eyewitness observations, we apply the ‘Tsunami Squares’ method to model the Hongyanzi landslide and its generated waves. Landslide simulations indicate a maximum impact velocity of ∼16 m/s that matches well with an eyewitness video. The computed post-slide mass stopped on the near riverbed with a shape fitting the observed geological profile. Tsunami simulations reveal a large region of wave impacts that coincide with the observed runup heights. The successful reproduction of the dynamics of this landslide-generated river tsunami emphasizes the capacity and efficiency of Tsunami Squares modeling in emergency reaction and risk assessment.〈/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: 〈span〉〈div〉Summary〈/div〉Improvement of global 3D Earth density and velocity models is based in part on measurements of Earth’s normal mode eigenfrequencies and splitting function coefficients. Despite many methods developed inconsistency in measurements still exists and it is difficult to understand which results are more precise, that is which methods introduce less systematic biases in the measurements. Therefore, the main goal of this study is to test the performances of typically used techniques in low-frequency normal mode studies: the optimal sequence estimation stacking method and the autoregressive method in the frequency domain, where validation of the estimates is performed with the phasor walkout method. Motivations for their utilization are their easy and fast implementation and their accurate performances when it comes to eigenfrequency estimates. For this purpose, we first perform the analysis with synthetic seismograms in order to evaluate how the station distributions and noise levels impact the estimates of eigenfrequencies and structure coefficients. Synthetic seismograms are calculated for a 3D realistic Earth model, which includes Earth’s rotation as well as ellipticity and other lateral heterogeneities. They were computed by means of normal mode summation and a perturbation theory for modes up to 1 mHz. The three methods above are also applied to long-period seismometer and superconducting gravimeter data recorded after six earthquakes of magnitude greater than 8.3. Finally, our study shows that the optimal sequence estimation is sensitive to the station distribution under the noise influence, while the autoregressive method for frequency estimation gives us reasonable estimates within the estimated error bars. Moreover, we present new estimates of eigenfrequencies and Q-factors for 〈sub〉0〈/sub〉S〈sub〉2, 0〈/sub〉S〈sub〉3, 2〈/sub〉S〈sub〉1〈/sub〉 and 〈sub〉3〈/sub〉S〈sub〉1〈/sub〉 multiplets. A new value for the 〈span〉c〈/span〉〈sub〉20〈/sub〉 structure coefficient of 〈sub〉0〈/sub〉S〈sub〉2〈/sub〉 multiplet −0.7233 ± 0.0623 μHz is obtained.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We propose a novel approach to compute the gravity field due to density anomaly in both outside and inside of the solid Earth with high accuracy and efficiency. The high accuracy comes from the direct employment of the analytic gravitation solution between any point on a two-dimensional (2D) plane in the horizontal direction and individual mass cubes. The high computational efficiency comes from two aspects: 1, the application of the highly efficient 2D discrete convolution algorithm; and 2, a newly developed algorithm for the optimized computation of the weight coefficient matrix. Numerical examples for applying to compare with analytical solutions demonstrated its excellent accuracy. Comparison with other state-of-the-art gravity modeling algorithm has proved that this algorithm has superior performance in both accuracy and efficiency. Application to analyze real topography demonstrated the practicality. This algorithm will be an attractive candidate for carrying out the forward modeling step in geophysical inversion problems with the claimed and proved advantages.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We assess the stress field that drives lithospheric deformation, focusing on the active New Zealand plate boundary zone between the Australian and Pacific plates. Here there is a rich database for the horizontal velocity field and crustal structure, used to derive gravitational potential energy (GPE). We solve the stress balance equations, in the context of a thin sheet model of a viscously deforming lithosphere, for characteristic deviatoric stresses and viscosities, defined as the integral of these with depth (divided by the layer thickness), using the stress method of Flesch 〈span〉et al.〈/span〉 (2001), where the input parameters are the fields of strain rate and GPE that apply to the sheet. Synthetic tests show that the stress method is able to resolve the stress field to high (20 per cent) levels of noise in the input strain rates, and the mean stress is a very robust feature of the inversions, regardless of noise levels. We invert for the stress and viscosity fields in New Zealand, calculating the field of GPE from topography/bathymetry and crustal data, and the field of strain rates from either a long term (multi millennial) velocity field inferred from the rate and pattern of Quaternary faulting, or a short term (decadal) velocity field directly observed with decadal Global Positioning System (GPS) measurements. In addition, we consider the effect of shear stresses on the subducted plate interface along the Hikurangi Margin (5–15 MPa), or regionally. We explore the effect of GPE on the inversion results by calculating these for a range of deforming layers (ie thin sheet with 35–150 km thickness), in effect sampling the lithospheric strength in the crust and mantle. The results show that the derived stress magnitudes (square root second invariant of the stress deviator) are in the range 0–35 MPa, with mean values of 13 ± 1 MPa for all models, comparable to typical earthquake stress drops. Gravitationally induced stresses account for approximately half of the full deviatoric stress. Effective characteristic viscosities are 0.5 - 5 × 10〈sup〉21〈/sup〉 Pas in the deforming zone, with an approximate inverse relation between strain rate and viscosity, most likely controlled by thermal structure and/or lithology.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉I present a source-independent fracture imaging method to use passive seismic data for mapping subwavelength natural fractures. Unlike conventional source-dependent imaging that often adopts reflection-type seismic imaging with known source that is not available in passive seismic surveys, the proposed fracture imaging approach relies on the transmission and diffraction data without the need for source information. I assume that passive seismic data can be decomposed into two types of data: primary transmission wave data and diffraction (coda) wave data. The imaging formula states that primary waves should coincide with coda waves at scatterer points at the time of scattering. Instead of generating source wavefields in the conventional imaging method, the proposed method only need to propagate transmission wave data and diffraction wave data from the receiver arrays and apply an imaging condition to produce an image of fractures. This imaging procedure can be used for processing P-wave or S-wave. In synthetic examples, I evaluate the proposed method in several aspects: inaccurate source location, inaccurate velocity model, sparse receivers and irregular receiver spacing, elastic data, and joint surface and borehole acquisitions. I found that the proposed approach performed well (or even better) comparable to source-dependent fracture imaging when assuming exact source information is known. With perturbed source locations with random shifts (e.g., estimated source location with errors), however, fractures were missing in the source-dependent fracture imaging results but the proposed approach was not influenced. In the presence of velocity errors and sparse and irregular receiver spacing, the proposed method produces better fracture images than the source-dependent imaging results.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We introduce a formalism for estimating local spatial averages of the core-mantle boundary (CMB) radial magnetic field and its time derivatives, based on magnetic field observations collected by low-Earth-orbit satellites. This provides a useful alternative to conventional core field modelling based on global spherical harmonic basis functions, where noise in the polar regions maps into all harmonics, and model regularization and spectral truncation are required. A powerful perspective offered by the proposed technique is formal appraisal of the spatial resolution and variance of the resulting field averages. We use the Green’s functions for the Neumann boundary value problem to link the satellite observations to the radial magnetic field on the CMB and estimate field averages using a modified Backus-Gilbert inversion approach. Our approach builds on the Subtractive Optimally Localized Averages (SOLA) method developed in helioseismology, that seeks averaging kernels as close as possible to a chosen target kernel. We are able to account for both internal and external field sources and can easily incorporate data error covariance information, for example describing along-track serial error correlation. As a proof of concept we present a global map collecting local estimates of the radial main field (MF) constructed on a grid at the CMB with one degree spacing in latitude and longitude, derived from one month of three component vector magnetic field data collected by the 〈span〉Swarm〈/span〉 satellite trio, using data from dark and geomagnetically quiet times. Using sums and differences of the field components taken along track and in the east-west direction we obtain estimates with spatial resolution kernel widths varying between 18 and 54 degrees depending on the latitude, and a standard deviation of approximately 10μT (i.e. 5% of the mean CMB field amplitude). The morphology of our CMB radial field map agrees well with results from conventional spherical harmonic field models. In a second application, we determine local estimates of the average rate of change, or secular variation (SV), of the radial field at the CMB, initially considering two year time windows, and performing the analysis on data collected by either the 〈span〉Swarm〈/span〉 or CHAMP satellites. We obtain stable local estimates of the SV at the CMB, and present maps of estimates with averaging kernel widths of approximately 42, 33 and 30 degrees on the equator, with corresponding standard derivations of 0.25, 2.5 and 5 μT/yr. By subtracting SV estimates constructed at different epochs we are able to calculate the local aggregated secular acceleration (SA) and to study its time changes. Differencing SV estimates 2 years apart, and considering an averaging kernel width of 42 degrees on the equator, we obtain SA maps very similar to those found in the CHAOS-6-x7 field model truncated at SH degree 10. Using our approach we are able to directly control the width of the spatial averaging kernel and the length of the time window, enabling us to directly study the robustness of the inferred SA. Pushing to higher resolution in time, considering one year differences of SV estimates constructed using one year windows, we are able to track the evolution of coherent SA structures in time-longitude plots at the equator. At 25° W in mid 2007 we find a distinctive SA ’cross-over’ event, with strong, oppositely signed and adjacent, SA features rapidly changing sign within a year. Our method is well suited for studying such spatio-temporally localized SA events at high resolution; there will be further opportunities for such investigations as the time series of data provided by the 〈span〉Swarm〈/span〉 mission lengthens.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The evolution of the ratio between P- and S-wave velocity (〈span〉V〈/span〉〈sub〉P〈/sub〉/〈span〉V〈/span〉〈sub〉S〈/sub〉) with increasing fluid-saturated porosity is computed for isotropic rocks containing spheroidal pores. The ratio 〈span〉V〈/span〉〈sub〉P〈/sub〉/〈span〉V〈/span〉〈sub〉S〈/sub〉 is shown to either decrease or increase with increasing porosity, depending on the aspect ratio α of the pores, fluid to solid bulk modulus ratio ζ, and Poisson’s ratio ν〈sub〉0〈/sub〉 of the solid constituents of the rock. A critical initial Poisson’s ratio ν〈sub〉0, crit〈/sub〉 is computed, separating cases where 〈span〉V〈/span〉〈sub〉P〈/sub〉/〈span〉V〈/span〉〈sub〉S〈/sub〉 increases (if ν〈sub〉0〈/sub〉 〈 ν〈sub〉0, crit〈/sub〉) or 〈span〉decreases〈/span〉 (if ν〈sub〉0〈/sub〉 〉 ν〈sub〉0, crit〈/sub〉) with increasing porosity. For thin cracks and highly compressible fluids, ν〈sub〉0, crit〈/sub〉 is approximated by 0.157 ζ/α, whereas for spherical pores ν〈sub〉0, crit〈/sub〉 is given by 0.2 + 0.8ζ. When ν〈sub〉0〈/sub〉 is close to ν〈sub〉0, crit〈/sub〉, the evolution of 〈span〉V〈/span〉〈sub〉P〈/sub〉/〈span〉V〈/span〉〈sub〉S〈/sub〉 with increasing fluid-saturated porosity is near neutral and depends on subtle changes in pore shape and fluid properties. This regime is found to be relevant to partially dehydrated serpentinites in subduction zones (porosity of aspect ratio near 0.1 and ζ in the range 0.01–0.1), and makes detection of these rocks and possibly elevated fluid pressures difficult from 〈span〉V〈/span〉〈sub〉P〈/sub〉/〈span〉V〈/span〉〈sub〉S〈/sub〉 only.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Secondary microseismic sources emit seismic waves over long time spans. Reoccurring signals with similar slowness and frequency therefore arrive at seismic arrays. Blind source separation techniques can be used to identify and isolate such reoccurring signals from other signals and from diffuse seismic noise. Along these lines, we use non-negative matrix factorization as blind source separation technique to decompose continuous seismic array records. We model the recorded energy as a mixture of a few components with static slowness-frequency and time dependent amplitudes. Components and amplitudes are fitted to optimally explain the recorded seismic energy over time. These components represent secondary microseismic signals with quasi-static slowness-frequency vector and fluctuating amplitude. Each fitted component reveals the geographical origin (through the slowness-frequency vector) and time evolution of an active secondary microseism with high precision because it is separated from other signals and diffuse seismic noise. Furthermore, relative travel times can be automatically extracted for the signals that correspond to a specific component that can potentially be used in tomographic studies. We show two examples of seismic signals that were extracted with this technique, one focusing on P-waves from the typhoons Goni and Atsani, and another showing secondary microseism PKP signals from typhoon Glenda.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Locating microseismic events is essential for many areas of seismology including volcano and earthquake monitoring and reservoir engineering. Due to the large number of microseismic events in these settings, an automated seismic location method is required to perform real time seismic monitoring. The measurement environment requires a precise and noise-resistant event location method for seismic monitoring. In this paper, we apply Multichannel Coherency Migration (MCM) to automatically locate microseismic events of induced and volcano-tectonic seismicity using sparse and irregular monitoring arrays. Compared to other migration-based methods, in spite of the often sparse and irregular distribution of the monitoring arrays, the MCM can show better location performance and obtain more consistent location results with the catalogue obtained by manual picking. Our MCM method successfully locates many triggered volcano-tectonic events with local magnitude smaller that 0, which demonstrates its applicability on locating very small earthquakes. Our synthetic event location example at a carbon capture and storage site shows that continuous and coherent drilling noise in industrial settings will pose great challenges for source imaging. However, automatic quality control techniques including filtering in the frequency domain and weighting are used to automatically select high quality data, and can thus effectively reduce the effects of continuous drilling noise and improve source imaging quality. The location performance of the MCM method for synthetic and real microseismic datasets demonstrates that the MCM method can perform as a reliable and automatic seismic waveform analysis tool to locate microseismic events.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The unification of local height systems has been a classical geodetic problem for a long time, the main challenges of which are the estimation of offsets between different height systems and the correction of tilts along the levelling lines. It has been proposed to address these challenges with clock networks. The latest generation of optical clocks as well as the dedicated frequency links, e.g., optical fibres, are now approaching to deliver the comparison of frequencies at the level of 1.0 × 10〈sup〉−18〈/sup〉. It corresponds to an accuracy of about 1.0 cm in height difference. Clock networks can thus serve as a powerful tool to connect local height systems. To verify the idea, we carried out simulations using the EUVN/2000 (European Unified Vertical Network) as a priori input. Four local height systems were simulated from the EUVN/2000 by introducing individual offsets and tilts, and were re-unified by using measurements in clock networks. The results demonstrate the great potential of clock networks for height system unification. In case that the offsets between different height systems and tilts along national levelling lines in both longitudinal and latitudinal directions are considered, three or four clocks measurements for each local region are sufficient for the unification. These clocks are to be interconnected and should be properly arranged so that they can sense the levelling tilts where necessary. Our results also indicate that even clocks with one magnitude poorer accuracy than the desired ones can still unify the height systems to some extent, but it may cause a shift for the re-unified system.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We investigate slip-distribution models of the 2011 Tohoku earthquake, with a particular focus on diffracted tsunamis and uplift-induced waves along the back-arc region of the Japanese Island Arc.The 2011 Tohoku earthquake produced a large amplitude tsunami that diffracted around Kyushu Island before reaching Korea. At the same time, this earthquake co-seismically induced short-period small-amplitude sea waves in the East Sea. We performed tsunami simulations using seven fault models of the Tohoku earthquake to examine whether the models can accurately reproduce the observed waveforms in the open sea of the western Pacific Ocean, the South Sea of Korea, and the coast of the East Sea. For each fault model, we investigate tsunami features due to geomorphological characteristics of the Korean Peninsula in the Korea offshore. To determine which slip distribution model shows a good performance in the tsunami simulations, we set three criteria; the delay time between observations and synthetic waveforms, the normalized mean residual, and the normalized RMS misfit. Depending on the study region, all models show varying degrees of accuracy. The fault dimensions and the amount of slip have a larger effect on the RMS misfit then the slip distribution patterns of the fault models for observations along the Korean coast and the western coast of Japan.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The study of radiation patterns in poroelastic media allows us to visually explore the possibility of reconstructing crucial model properties describing reservoir rocks, and to examine the coupling effects between different parameters during full-waveform inversion (FWI). In this paper, we derive analytical formulae for the radiation patterns of single parameter perturbations in fluid-saturated porous media by deriving scattered wavefields based on plane-wave theory and the far-field approximation. We illustrate these scattered wavefields via their radiation patterns expressed as a function of the angle between the incident and scattered waves. To simplify the algebra, we consider poroelastic waves at seismic (low) frequencies, where the fast compressional wave and shear wave are propagating modes but the slow compressional wave is severely dispersive. To verify our derivation of the analytical radiation patterns, we also compute them numerically by perturbing one parameter at a single point, keeping the other parameters fixed at their background values. We find that all the analytical radiation patterns match the wavefronts of the numerically computed scattered wavefields, well indicating that our derivations are correct. Parameters such as the solid density, fluid density, viscosity of the fluid, and intrinsic permeability, have similar radiation patterns and thus show strong coupling effects. Therefore, we anticipate difficulties in recovering these parameters in a multi-parameter FWI procedure. In an attempt to mitigate these trade-offs, we analyze different parameterizations which result in different radiation patterns. As the patterns that we observe are similar to those of the elastic case, we anticipate that parameter separation might be easier when inverting for velocities than for moduli.〈/span〉

Description: 〈span〉〈div〉Abstract〈/div〉In the Sea of Marmara, areas of gas seepage or cold seeps are tightly related to the faults system and understanding the spatial and temporal dynamics in gas-related processes is crucial for geohazard mitigation. Although acoustic surveys proved to be efficient in detecting and locating cold seeps, temporal variability or trends in the gas-related processes are still poorly understood. Two arrays of 10 ocean bottom seismometers were deployed in the western part of the Sea of Marmara in 2011 and 2014, respectively. In addition to the local seismic events, the instruments recorded a large number of short duration events and long-lasting tremors. Short duration events are impulsive signals with duration 〈 1 s, amplitude well above the noise level and a frequency spectrum with one or two narrow peaks. They are not correlated from one site to another, suggesting a very local source. Tremors consist of sequences of clustered impulsive signals lasting for minutes to more than an hour with a multi-peak frequency spectrum. Based on evidence of known seepage and by analogy with volcanic and hydrothermal models, we suggest that short duration events and tremors are associated with gas migration and seepage. There is a relationship between tremors associated with gas emission and the local seismicity, although not systematic. Rather than triggering gas migration out of the seabed, locally strong earthquakes act as catalysts when gas is already present or gas emission is already initiated.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Long-period Rayleigh wave Horizontal to Vertical amplitude (H/V) ratios at a station provide information about local earth structure that is complementary to phase velocity. However, a number of studies (Ferreira and Woodhouse, 2007; Lin 〈span〉et al.〈/span〉, 2012) have observed that significant scatter appears in these measurements making it difficult to use H/V ratio measurements to resolve earth structure. Some of the scatter in these measurements has been attributed to local geological structure while some has remained unaccounted for. Most Global Seismographic Network (GSN) stations contain two nearby high-quality broadband seismometers (e.g. in the same vault, but on different piers or in different boreholes). For each broadband sensor in the IRIS/USGS component of the GSN, we estimate H/V ratios of fundamental mode Rayleigh waves using M 〉 6.5 earthquakes from 2001 to 2018 (around 19,000 measurements). We compute these ratios at a number of discrete periods (25, 50, 75, 100, and 150 s) and find that for well-isolated Rayleigh waves (windows where the correlation coefficients between radial and the phase-shifted vertical components are greater than 0.9) significant scatter in H/V ratios occurs between co-located sensors (greater than 25 per cent at 100 s period). This suggests the scatter in H/V ratio measurements can be at least partially attributed to extremely local phenomena such as sensor emplacement in the vault. We also find that H/V ratios can vary as a function of event back-azimuth, indicating that care must be taken when computing average ratios for a station, as a large number of events from a given region could bias H/V ratio measurements at a station.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We analyse noise characteristics of spatial and temporal correlation of 260 continuous GPS sites from Crustal Movement Observation Network of China (CMONOC). These data sets were mainly collected between 2010 and 2016, with an average of 6 yr of position time-series. In the functional analysis, a clear regional dependence of seasonal movements has been observed and other significant periodic signals are detected nearby the GPS draconitic period and its harmonics. The distribution of these periodic signals shows a spatial correlation, along with non-negligible local inconsistencies. In the stochastic analysis, impacts of the periodicities on the noise assessment have been investigated and Maximum likelihood estimation is used to study noise properties for deseasonalized residual time-series having the seasonal signals removed and filtered residual time-series having other periodic signals removed further. We demonstrate that for both solutions, the flicker noise is thought to be the dominant time-correlated noise and velocity uncertainties may be underestimated 8–10 times if assuming a pure white noise. Ignoring the periodicities could bias the estimation of noise amplitude, spectral index and velocity uncertainty. After removing the periodic signals, the median flicker noise magnitude shows an average of ∼10 per cent reductions and the noise process shifts closer towards white noise. A correlation between the index variations and RMS variations has been observed, indicating that the index varies more significantly for sites with more periodic signals removed. Besides, the spectral index in the vertical component has a better spatial correlation than that in the horizontal and the spatial distribution of the index of deseasonalized solutions seems to correlate well with the amplitudes of seasonal signals, probably implying common sources of spatial variations of these characteristics. Furthermore, analyses of intersite correlations indicate that the correlation induced from the periodic signals displays a similar pattern to the deseasonalized solutions, confirming that the period pattern is spatially correlated and can induce time-series correlations. Finally, the stochastic processes of the common mode noise are predominantly featured by spatially correlated flicker and white noise over a wide range, consistent with previous results.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉As temperature increases with depth and the creep resistance of rock decreases exponentially, a high-viscosity sub-lithospheric layer, just beneath the ‘elastic’ lithosphere is expected to exist. Depending on the temperature profile, a low-viscosity asthenosphere may also exist if the temperature deeper down gets high enough. Since the temperature profile is expected to change laterally – especially from below the oceans to cratonic areas underneath continents, rock properties of the lithosphere, high-viscosity sub-lithosphere and low-viscosity asthenosphere are expected to change laterally. Our aim is to constrain sub-lithospheric properties (depth, thickness and viscosity), lateral lithospheric thickness variations and asthenospheric properties using observed GIA data. A Coupled Laplace-Finite Element Method is used to compute gravitationally self-consistent sea level with time-dependent coastline and rotational feedback in addition to changes in deformation, gravity and the state of stress. We start with the VM5a-ICE-6G_C model combination and then modify the lithospheric, sub-lithospheric and asthenospheric properties (including lateral thickness variation) while keeping the mantle viscosities the same as VM5a. Through this study, we confirm that the sub-lithospheric and asthenospheric properties can significantly affect the predicted global relative sea level (RSL), present-day gravity rate-of-change (g-dot) and uplift rate (u-dot) in Laurentia and Fennoscandia. In addition, incorporating the elastic lithosphere with lateral thickness variation, sub-lithosphere and asthenosphere can improve the fit to global RSL, but the predicted peak values of g-dot and u-dot in Laurentia may decrease slightly but not significant enough to affect the fit to the observed data. Our results prefer an elastic lithosphere that has maximum thickness of 140 km under continental cratons but reduces to 60 km underneath the oceans. The results preferred depth of the asthenospheric bottom is around 190–200 km with asthenospheric viscosity around 10〈sup〉20〈/sup〉Pa s. Finally, we show that the best laterally heterogeneous mantle model we found in previous publication when combined with the lithosphere with lateral thickness variaion gives the best fit to global RSL and peak g-dot and u-dot in Laurentia simultaneously.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We implement an analytical model based on flexural deflection of a thin elastic disc to investigate the magnitude of lithospheric decompression caused by deglaciation at upper crustal magmatic reservoirs. Considering a published numerical climate model describing the space–time evolution of deglaciation after the Last Glacial Maximum (LGM) along the Southern Volcanic Zone (SVZ) of the Andes, we demonstrate that changes in pressure at upper crustal levels (〈10 km depth) at the scale of several hundred years are of the order of 10–100 MPa. Total decompression and decompression rate (300–150 kPa yr〈sup〉−1〈/sup〉) are 1–2 orders of magnitude larger than values previously estimated by other authors who assume that glacial loads are supported by an elastic half-space, that is, of infinite elastic thickness. The large decompression caused by flexural unbending of an elastic plate of finite thickness as assumed here can easily surpass the tensile strength of rocks (5–20 MPa), creating adequate conditions for failure of the reservoir walls, dike propagation inside and outside the reservoir and the eventual collapse of the reservoir accompanying an explosive eruption. We apply our results to the analysis of post-glacial eruptions of SVZ volcanoes, which erupted large volumes (〉10 km〈sup〉3〈/sup〉) of mafic ignimbrites hundreds to thousands of years after deglaciation onset. We show that this time lag is necessary to achieve a decompression of several tens of megapascals at depths of several kilometres that are consistent with the location of magmatic reservoirs as estimated by independent petrologic, seismic and/or geodetic studies. Moreover the northward increase of this time lag is in agreement with a smaller size of the Andean ice cap in the north than in the south during the LGM. For wet, volatile-rich magmas typical of subduction zones, the effect of large decompression at upper crustal reservoirs caused by flexural unbending of the lithosphere after deglaciation could play a major role in promoting large explosive eruptions through devolatization of the magma, during past deglaciation events as demonstrated here for the LGM along the SVZ and current accelerated ice retreat caused by climate change over large segments of subduction-related arcs at higher latitudes.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The unification of local height systems has been a classical geodetic problem for a long time, the main challenges of which are the estimation of offsets between different height systems and the correction of tilts along the levelling lines. It has been proposed to address these challenges with clock networks. The latest generation of optical clocks as well as the dedicated frequency links, for example optical fibres, are now approaching to deliver the comparison of frequencies at the level of 1.0 × 10〈sup〉−18〈/sup〉. It corresponds to an accuracy of about 1.0 cm in height difference. Clock networks can thus serve as a powerful tool to connect local height systems. To verify the idea, we carried out simulations using the EUVN/2000 (European Unified Vertical Network) as 〈span〉apriori〈/span〉 input. Four local height systems were simulated from the EUVN/2000 by introducing individual offsets and tilts, and were reunified by using measurements in clock networks. The results demonstrate the great potential of clock networks for height system unification. In case that the offsets between different height systems and tilts along national levelling lines in both longitudinal and latitudinal directions are considered, three or four clocks measurements for each local region are sufficient for the unification. These clocks are to be interconnected and should be properly arranged so that they can sense the levelling tilts where necessary. Our results also indicate that even clocks with one magnitude poorer accuracy than the desired ones can still unify the height systems to some extent, but it may cause a shift for the reunified system.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Waveform backprojection (BP) is a key technique of earthquake-source imaging, which has been widely used for extracting information of earthquake source evolution that cannot be obtained by kinematic source inversion. The technique enjoys considerable popularity, owing to the simplicity of its implementation and the robustness of its processing, but the physical meaning of BP images has remained elusive. In this study, we reviewed the mathematical representation of BP and hybrid BP (HBP) methods, following the pioneering work of Fukahata 〈span〉et al〈/span〉. (〈a href="http://academic.oup.com/gji#bib15"〉2014〈/a〉), to clarify the physical implications of BP images. We found that signal intensity in BP and HBP images is scaled with the amplitude of the Green’s function that corresponds to a unit-step slip, which results in the signal intensity being depth dependent. We propose variants of BP and HBP, which we call kinematic BP and HBP, respectively, to relate the BP signal intensity to slip motion of an earthquake by modifying the normalizing factors used in the original BP and HBP methods. The original BP and HBP images remain useful for assessing the spatiotemporal strength of the wave radiation, which scales with the amplitude of the Green’s function, whereas the kinematic BP and HBP methods are suitable for imaging the slip motion that is responsible for the high-frequency radiation produced during the source-rupture process.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉 〈strong〉Moment Tensors of hydraulically induced AEs:〈/strong〉 Hydraulic fracturing is an important technique in the development of enhanced geothermal systems and unconventional resources. Although the fracture modes induced by hydraulic fracturing influence the recovery efficiency of the resources, the current understanding of this relationship is insufficient. In this study, we considered the acoustic emissions (AEs) induced during hydraulic fracturing under uniaxial loading conditions in the laboratory, and applied a moment tensor analysis by carefully correcting the coupling condition and directivity of AE transducers. Experiments were conducted for two types of Kurokami–jima granite samples: those with a rift plane perpendicular (Type H) or parallel (Type V) to the expected direction of fracture propagation (i.e. along the loading axis). In the experiments, both sample types experienced a significant number of shear, tensile and compressive events. The dominant fracture mode for Type H samples is found to be tensile events in which the fracture plane is parallel to the loading axis, whereas for Type V samples, shear events are dominant. This difference suggests that the dominant fracture modes induced by hydraulic fracturing are highly dependent on the relationship between the direction of fracture propagation and orientation of pre-existing weak planes.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We present a sensitivity analysis aimed at testing whether observables related to glacial isostatic adjustment can support or refute the occurrence of a low viscosity melt-rich layer (MRL) above the mantle transition zone, as required by the ‘transition-zone water-filter’ model (Bercovici & Karato 2003). In total, 1600 model runs were performed sampling a range of MRL thicknesses (1, 10 and 20 km) and viscosities (10〈sup〉15〈/sup〉–10〈sup〉19〈/sup〉 Pa·s), plausible viscosity values in the upper and lower mantle regions and four distinct ice histories. To determine decay time constraints, we consider relative sea level (RSL) data from two sites [Ångerman River (ÅR), Sweden and Richmond Gulf (RG), Canada] and use a new method of observational sea level data correction. Comparing model output of postglacial decay times and ${\skew{6}\dot{J}_2}$ to observational constraints, we find numerous possible solutions, largely as a result of parameter trade-off. The investigated observables are sensitive to the existence of an MRL and reasonable variations in its thickness and viscosity. The magnitude and nature of this sensitivity varies between the two data types as well as the adopted background viscosity structure. Decay time results from either considered location do not strictly support or exclude MRL existence. However, both locations offer MRL viscosity requirements for given thicknesses, with ÅR being more restrictive. RG constraints allow MRL viscosities as low as 10〈sup〉16〈/sup〉 Pa·s (10 km) and 10〈sup〉17〈/sup〉 Pa·s (20 km). ÅR results narrow these permitted viscosity ranges to 10〈sup〉18〈/sup〉 Pa·s or greater for both 10 and 20 km MRL thicknesses. In the case of a 1 km thick MRL, ÅR constraints permit the viscosity to be as low as 10〈sup〉17〈/sup〉 Pa·s, whereas those of RG permit any MRL viscosity. The decay time observations are satisfied by only a small subset of ‘background’ mantle viscosities (regardless of the MRL properties), none of which support a spherically symmetric solution of Earth viscosity. Finally, comparing model output to the observed ${\skew{6}\dot{J}_2}$ value did not provide any constraints on MRL properties. However, our results show that this observable has a strong preference for viscosity values in the lower mantle that are equal to or greater than 10〈sup〉22〈/sup〉 Pa·s.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We investigated an 〈span〉M〈/span〉〈sub〉w〈/sub〉 ∼ 6.2 earthquake doublet on the border of the USA and Canada using ALOS2 Light-of-Sight displacements and GPS measurements. We selected three L-band ALOS-2 interfergorams with temporal baselines of one yr to extract coseismic deformation maps, in which master and slave images were both acquired in July. A subpixel-based alignment and another range spectral splitting techniques under the GAMMA InSAR software framework were applied to improve the interferometric coherence and reduce the effects of phase anomalies in two of the three interferometric pairs due to either ionospheric delay or a potential focusing issues in the generation of the ALOS2 SLC data. The updated interferograms convincingly reveal deformation fringe patterns produced by the two earthquakes. We conducted a nonlinear geophysical inversion to estimate the geometric parameters of the earthquakes with the InSAR and GPS measurements. The best-fitting model shows that a thrust faulting on a reverse fault and left-lateral strike-slip faulting on a nearly vertical fault with the centroid depths of 9.3±0.6 and 8.4±0.7 km, respectively, are most likely responsible for the earthquake doublet. The eastern Denali fault (EDF) and Duke River fault are major active faults in the region and the earthquake doublet could be due to reactivation of the part of the two faults system.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉In this paper, we present a series of mathematical abstractions for seismologically relevant wave equations discretized using finite-element methods, and demonstrate how these abstractions can be implemented efficiently in computer code. Our motivation is to mitigate the combinatorial complexity present when considering geophysical waveform modelling and inversion, where a variety of spatial discretizations, material models, and boundary conditions must be considered simultaneously. We accomplish this goal by first considering three distinct classes of abstract mathematical models: (1) those representing the physics of an underlying wave equation, (2) those describing the discretization of the chosen equation onto a finite-dimensional basis and (3) those describing any spatial transforms. A full representation of the discrete wave equation can then be constructed using a hierarchical nesting of models from each class. Additionally, each class is functionally orthogonal to the others, and with certain restrictions models within one class can be interchanged independently from changes in another. We then show how this recasting of the relevant equations can be implemented concisely in computer software using an abstract object-oriented design, and discuss how recent developments in the numerical and computational sciences can be naturally incorporated. This builds to a set of results where we demonstrate how the developments presented can lead to an implementation capable of multiphysics waveform simulations in completely unstructured domains, on both hypercubical and simplical spectral-element meshes, in both two and three dimensions, while remaining concise, efficient and maintainable.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Traveltime approximation plays an important role in seismic data processing, for example, anisotropic parameter estimation and seismic imaging. By exploiting seismic traveltimes, it is possible to improve the accuracy of anisotropic parameter estimation and the resolution of seismic imaging. Conventionally, the traveltime approximations in anisotropic media are obtained by expanding the anisotropic eikonal equation in terms of the anisotropic parameters and the elliptically anisotropic eikonal equation based on perturbation theory. Such an expansion assumes a small perturbation and weak anisotropy. In a realistic medium, however, the assumption of small perturbation likely breaks down. We present a retrieved zero-order deformation equation that creates a map from the anisotropic eikonal equation to a linearized partial differential equation system based on the homotopy analysis method. By choosing the linear and nonlinear operators in the retrieved zero-order deformation equation, we develop new traveltime approximations that allow us to compute the traveltimes for a medium of arbitrarily strength anisotropy. A comparison of the traveltimes and their errors from the homotopy analysis method and from the perturbation method suggests that the traveltime approximations provide a more reliable result in strongly anisotropic media.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Guided waves in a water layer overlaying an elastic half-space are known as normal modes. They are often present in seismic recordings at long offsets in shallow-water environment and generally considered coherent noise. The normal modes, however, carry important information about the near-surface and, as demonstrated by a number of authors, can be used to obtain the shallow velocity model. There is a growing evidence that the latter needs not to be isotropic due to various geological reasons. Motivated by that, we consider the normal-mode propagation in case the elastic half-space exhibits orthorhombic anisotropy. We derive the period equation that describes the normal-mode phase velocity dispersion. To simplify the complicated expression, we present acoustic and ellipsoidal orthorhombic approximations. We also outline the approach towards the group velocity and group azimuth calculation and apply it to the ellipsoidal case to obtain concise and intuitive expressions. Using numerical test, we study the relation between phase and group domains in elastic orthorhombic case. The deviation between velocities and azimuths in these domains is the strongest for low frequencies and it rapidly decreases with increasing frequency. For higher frequencies, the anisotropy effects of the underlaying half-space are barely detectable since the observed signal is composed mainly of the direct acoustic wave, resulting in the two domains being nearly indistinguishable.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Understanding why Earth’s lithosphere is divided into several plates while other terrestrial bodies have unbroken lids is a long-standing challenge, often addressed with the help of numerical modelling. A key mechanism defining the transition between these two convective regimes is the formation of shear zones that cut through the entire lithosphere in regions with high stresses. Here we present a modelling study in which lithospheric stresses resulting from small-scale convection in the upper mantle are analysed. We perform model simulations that include elasticity and a free surface and evaluate how these physical complexities affect stress distribution inside the lithosphere, which in turn controls the depths of yielding and the possible initiation of subduction. We show that the spatial distribution of stress is significantly altered by the presence of elastic deformation only when the model lithosphere acts as a thick plate capable of bending. Whether or not this is the case depends on the viscosity model. For an Arrhenius viscosity limited by a cut-off value that produces an essentially rigid lid, flexure dominates the observed lithospheric stress pattern in simulations with a free surface. The amplitudes of the stress are, when a free surface is assumed but elasticity is neglected, largely overestimated. Including both a free surface and elasticity results in stresses with maximum amplitudes close to those observed in the traditional models with a viscous rheology and a free-slip upper boundary, suggesting that having no additional complexity is, in a way, better than employing just a free surface. We also demonstrate how the use of impermeable free-slip side boundaries can result in the formation of unnatural, laterally locked convection cells, and bias the results of a parametric study. For each point in the parameter space, we perform several simulations with slightly different initial temperature fields in order to statistically eliminate the occurrence of locked states.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The design of an array configuration is an important task in array seismology during experiment planning. Often the array response function (ARF), which depends on the relative position of array stations and frequency content of the incoming signals, is used as the array design criterion. In practice, additional constraints and parameters have to be taken into account, for example, land ownership, site-specific noise levels or characteristics of the seismic sources under investigation. In this study, a flexible array design framework is introduced that implements a customizable scenario modelling and optimization scheme by making use of synthetic seismograms. Using synthetic seismograms to evaluate array performance makes it possible to consider additional constraints. We suggest to use synthetic array beamforming as an array design criterion instead of the ARF. The objective function of the optimization scheme is defined according to the monitoring goals, and may consist of a number of subfunctions. The array design framework is exemplified by designing a seven-station small-scale array to monitor earthquake swarm activity in Northwest Bohemia/Vogtland in central Europe. Two subfunctions are introduced to verify the accuracy of horizontal slowness estimation; one to suppress aliasing effects due to possible secondary lobes of synthetic array beamforming calculated in horizontal slowness space and the other to reduce the event’s mislocation caused by miscalculation of the horizontal slowness vector. Subsequently, a weighting technique is applied to combine the subfunctions into one single scalar objective function to use in the optimization process.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We develop a new method for measuring ellipticity of Rayleigh waves from ambient noise records by degree-of-polarization (DOP) analysis. The new method, named DOP-E, shows a good capability to retrieve accurate ellipticity curves separated from incoherent noise. In order to validate the method we perform synthetic tests simulating noise in a 1-D earth model. We also perform measurements on real data from Antarctica and Northern Italy. Observed curves show a good fit with measurements from earthquake records and with theoretical ellipticity curves. The inversion of real data measurements for 〈span〉vS〈/span〉 structure shows a good agreement with previous models. In particular, the shear-wave structure beneath Concordia station shows no evidence of a significant layer of liquid water at the base of the ice. The new method can be used to measure ellipticity at high frequency and therefore it will allow the imaging of near-surface structure, and possibly of temporal changes in subsurface properties. It promises to be useful to study near-surface processes in a wide range of geological settings, such as volcanoes, fault zones and glaciers.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The interaction of subducted oceanic lithosphere with the discontinuities of the mantle transition zone (MTZ) provides insight into the composition and temperature of the subducted slab as well as potential melting of the slab or the surrounding mantle and loss of volatiles from the slab. Detailed mapping of the structure of the MTZ will help to better understand how slabs transport material and volatiles into the mantle and how phase transitions affect the slab dynamics. Here we use a dense network of seismic stations in northern Anatolia to image the structure of the MTZ discontinuities in detail using 〈span〉P〈/span〉-wave receiver functions. With a station spacing of about 7 km and a surface footprint of ∼35 km × ∼70 km, analysing receiver functions calculated from teleseismic earthquakes that occurred during an ∼18-month deployment produced clear images of where the MTZ interacts with the Tethys/Cyprus slabs that either lie flat on the 660-km discontinuity or pass into the lower mantle. We observe an undulating 660-km discontinuity depressed by up to 30 km and a slightly depressed (1–2 km) 410-km discontinuity, apparently undisturbed by the slab. The MTZ is thickened to ∼270 km as result of the cool slab in the MTZ influencing the 660-km discontinuity and includes an arrival at ∼520-km depth likely from the top of a flat lying slab or a discontinuity related to a solid–solid phase transition in the olivine component of the mantle. We find evidence for low-velocity zones both above and below the 410-km discontinuity and above the 660-km discontinuity. The low-velocity zones around the 410-km discontinuity might be the result of hydration of the MTZ from the slab and upward convection of MTZ material into the upper mantle. The origin of the low-velocity zone around the 660-km discontinuity is less clear and could be related to sedimentation of subducted mid-ocean ridge basalts. The small footprint of the seismic array provides accurate information on the structure of the MTZ in an area influenced by subduction and shows small-scale changes in MTZ structure that might be lost in studies covering larger areas with sparser sampling.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We perform non-linear time-series analysis on a harmonic tremor seismogram recorded at 830 m away from the centre of the crater during the 2011 eruption at Shinmoedake, Japan. We found features suggesting the existence of period doubling bifurcation in the harmonic tremor signal, implying that the harmonic tremor might be generated by a non-linear process. In order to quantify the non-linearity in the harmonic tremor signal, we measure the correlation dimension 〈span〉D〈/span〉 and the maximal Lyapunov exponent λ. For one short but stable segment of the harmonic tremor seismogram, we obtained 〈span〉D〈/span〉 = 1.12 and λ = 0.03 s〈sup〉−1〈/sup〉. This result implies that the stable oscillation of the harmonic tremor is predominantly a limit cycle with small amounts of chaos present. We then use surrogate data analysis to check that our measurements of 〈span〉D〈/span〉 and λ do not include any false positive detection of non-linearity. Limit cycles imply that the harmonic tremor is generated by self-sustained oscillations. We show that the autonomous Julian tremor model is able to exhibit period doubling bifurcation. We also show that the non-autonomous Julian tremor model with a step like increase and decrease in input pressure is able to exhibit oscillations of varying amplitude while keeping a constant frequency spectrum. Both phenomena were observed at Shinmoedake. We also demonstrate that the non-autonomous Julian tremor model with a transient input pressure is able to exhibit long period-like events in addition to harmonic tremor-like events, implying that the same non-linear mechanism could be responsible for the generation of both type of events.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Reservoir models are numerical representations of the subsurface petrophysical properties such as porosity, volume of minerals and fluid saturations. These are often derived from elastic models inferred from seismic inversion in a two-step approach: first, seismic reflection data are inverted for the elastic properties of interest (such as density, 〈span〉P〈/span〉-wave and 〈span〉S〈/span〉-wave velocities); these are then used as constraining properties to model the subsurface petrophysical variables. The sequential approach does not ensure a proper propagation of uncertainty throughout the entire geo-modelling workflow as it does not describe a direct link between the observed seismic data and the resulting petrophysical models. Rock physics models link the two domains. We propose to integrate seismic and rock physics modelling into an iterative geostatistical seismic inversion methodology. The proposed method allows the direct inference of the porosity, volume of shale and fluid saturations by simultaneously integrating well-logs, seismic reflection data and rock physics model predictions. Stochastic sequential simulation is used as the perturbation technique of the model parameter space, a calibrated facies-dependent rock physics model links the elastic and the petrophysical domains and a global optimizer based on cross-over genetic algorithms ensures the convergence of the methodology from iteration to iteration. The method is applied to a 3-D volume extracted from a real reservoir data set of a North Sea reservoir and compared to a geostatistical seismic AVA.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We investigate slip-distribution models of the 2011 Tohoku earthquake, with a particular focus on diffracted tsunamis and uplift-induced waves along the backarc region of the Japanese Island Arc. The 2011 Tohoku earthquake produced a large amplitude tsunami that diffracted around Kyushu Island before reaching Korea. At the same time, this earthquake coseismically induced short-period small-amplitude sea waves in the East Sea. We performed tsunami simulations using seven fault models of the Tohoku earthquake to examine whether the models can accurately reproduce the observed waveforms in the open sea of the western Pacific Ocean, the South Sea of Korea, and the coast of the East Sea. For each fault model, we investigate tsunami features due to geomorphological characteristics of the Korean Peninsula in the Korea offshore. To determine which slip distribution model shows a good performance in the tsunami simulations, we set three criteria; the delay time between observations and synthetic waveforms, the normalized mean residual, and the normalized RMS misfit. Depending on the study region, all models show varying degrees of accuracy. The fault dimensions and the amount of slip have a larger effect on the RMS misfit then the slip distribution patterns of the fault models for observations along the Korean coast and the western coast of Japan.〈/span〉

Description: 〈span〉〈div〉Abstract〈/div〉Economically important deposits of nickel-copper and platinum group elements tend to be associated with small mafic-ultramafic intrusions characterized by compositional layering. Our rock magnetic study of an exploration drill core recovered from a typical small intrusion from the recently discovered Kun-Man'e ore field have found that concentration-independent hysteresis parameters of the studied rocks generally co-vary with petrological characteristics, such as major element abundances and compositions of olivine and clinopyroxene. Although monoclinic pyrrhotite and defect-poor magnetite both contribute to magnetic properties of the samples, we argue that hysteresis parameters are not affected significantly by pyrrhotite content. Rather, the variations in the hysteresis parameters reflect variations in median grain-size of magnetite crystals. The correlation between the rock-magnetic grain-size indicators and geochemical crystallization temperature indicators suggests that the size of magnetite crystals is controlled by cooling regime. Compared to conventional laborious crystal size distribution (CSD) studies, rock magnetic measurements can provide a quick and sensitive way to evaluate CSD of magmatic magnetite in mafic and ultramafic rocks.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The Reykjanes geothermal system is a high-temperature seawater system situated in SW-Iceland. Interferometric analysis of the Sentinel-1 satellite synthetic aperture radar (InSAR) data has been used to determine a time series of ground deformation induced by geothermal utilization between April 2015 and October 2017. Surface displacements have been estimated at coherent pixels, indicating a steady and linear subsidence within a sub-circular bowl centered on the well field at a maximum near-vertical rate of about 25 mm/yr, together with horizontal contraction. The average line-of-sight (LOS) displacement rates from ascending and descending tracks are inverted to determine the characteristics of the deformation source at depth, modeling the geothermal reservoir as a body of simple geometry within an elastic half space. The results indicate a deformation source at about 1 km depth contracting at a rate of (0.7–0.9) × 10〈sup〉5〈/sup〉 m〈sup〉3〈/sup〉/yr during the 2015–2017 period. Using pressure and temperature monitoring data at 900 m depth as well as an analysis of the reservoir structure and rock properties, we infer that the recent estimated volume change can be attributed to processes in the steam cap situated in the topmost part of the geothermal reservoir, in the 800–1200 m depth range. Processes involve a combination of compaction under pressure decrease and/or thermal contraction due to cooling of the rocks within or near the steam cap. The steam cap expanded in response to a sudden pressure drop resulting from the increase in extraction of geothermal fluids for a new power plant in 2006.〈/span〉

Description: 〈span〉〈div〉Abstract〈/div〉The seismic speed and anisotropy are two methods frequently used to image the assembly inside the Earth. We study the crust assembly beneath the Biga Peninsula and the surrounding area in Northwest Turkey using the accelerometer and broadband recordings where short-to-medium period (5-20 s) Love-Rayleigh surface waves are utilized to extract the group-phase speed data (fundamental mode). Single-station and two-station techniques are engaged to understand the detected surface waves for the speed and anisotropy assemblies. The single-station group speeds are inverted in a tomographic approach to attain the two-dimensional group speed diagrams. The least-squares inversion procedure is utilized to find the speed-depth profiles (one-dimensional) under each grid location. The one-dimensional results are cooperatively inferred to attain the three-dimensional appearance of the S-wave speeds below the measured region. This process is reiterated for Love and Rayleigh waves. Isotropic configuration is not sufficient to concurrently describe the present detected Love-Rayleigh surface waves. Vertical transverse isotropic crust assembly is found to better elucidate the detected data showing the Rayleigh-Love discrepancy. Complex arrangement of sills and dykes due to the widespread plutonic and volcanic activity in the region linked to the interaction between the Turkish plate and the African plate (northward subducting) is thought to depict the crust assembly deformations causing the detected long-wavelength vertical transverse isotropy. The mineral orientation within horizontal sills and vertical dykes following the magma flow, which is independent of seismic wavelength, adds to the detected anisotropy. The upper crust vertical transverse isotropy is mostly negative; i.e. SV-wave is faster than SH-wave, which is assumed to be due to the existence of dykes. The middle-to-lower crust vertical transverse isotropy is commonly positive; i.e. SH-wave is faster than SV-wave, which is assumed to be due to the existence of sills. The two-station analyses operating on cross-correlograms give analogous vertical transverse isotropic results to those of the single-station estimates.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The representative elementary volume (REV) is a critically important concept in fractured rock investigations as it tells us at what scale the fractured domain can be represented by an anisotropic tensor as opposed to requiring the details of each individual fracture for modeling purposes. Whereas the REV size and corresponding tensor characteristics for the hydraulic conductivity (〈span〉K〈/span〉) in fractured rock have been the subject of numerous previous investigations, no studies to date have focused on the electrical conductivity (σ). This is despite the fact that geoelectrical measurements are arguably the most popular means of geophysically investigating fractured rock, typically via azimuthal resistivity surveying where the observed electrical anisotropy is commonly used to infer hydraulic characteristics. In this paper, we attempt to fill this void and present a systematic numerical study of the impacts of changes in fracture-network properties on the REV size and equivalent tensor characteristics for both the electrical and hydraulic conductivities. We employ a combined statistical and numerical approach where the size of the REV is estimated from the conductivity variability observed across multiple stochastic fracture-network realizations for various domain sizes. Two important differences between fluid and electric current flow in fractured media are found to lead to significant differences in the REV size and tensor characteristics for σ and 〈span〉K〈/span〉; these are the greater importance of the matrix in the electrical case and the single-power instead of cubic dependence of electric current flow upon aperture. Specifically, the REV for the electrical conductivity will always be smaller than that for the hydraulic conductivity, and the corresponding equivalent tensor will exhibit less anisotropy, often with notably different principal orientations. These findings are of key importance for the eventual interpretation of geoelectrical measurements in fractured rock, where we conclude that extreme caution must be taken when attempting to make the link to hydraulic properties.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Stress is a crucial factor that affects the permeability, and the influence of stress on permeability has been studied extensively by many investigators. In this article the relationship between confining pressure and permeability is investigated. For this purpose, a void ratio- pressure model based on the theory of incompressible elastic porous media is first proposed, then, based on the model and Kozeny-Carman equation, a new pressure–permeability relationship is proposed. Finally the relationship is utilized to predict the permeability of coal and mudstone under different confining pressure, and these predictions are compared with other models and found to be in better agreement with experimental data reported recently.〈/span〉

Description: 〈span〉〈div〉Abstract〈/div〉We perform nonlinear time series analysis on a harmonic tremor seismogram recorded at 830 m away from the centre of the crater during the 2011 eruption at Shinmoedake, Japan. We found features suggesting the existence of period doubling bifurcation in the harmonic tremor signal, implying that the harmonic tremor might be generated by a nonlinear process. In order to quantify the nonlinearity in the harmonic tremor signal, we measure the correlation dimension 〈span〉D〈/span〉 and the maximal Lyapunov exponent λ. For one short but stable segment of the harmonic tremor seismogram, we obtained 〈span〉D〈/span〉 = 1.12 and λ = 0.03 s〈sup〉−1〈/sup〉. This result implies that the stable oscillation of the harmonic tremor is predominantly a limit cycle with small amounts of chaos present. We then use surrogate data analysis to check that our measurements of 〈span〉D〈/span〉 and λ do not include any false positive detection of nonlinearity. Limit cycles imply that the harmonic tremor is generated by self-sustained oscillations. We show that the autonomous Julian tremor model [〈span〉Julian〈/span〉, 1994] is able to exhibit period doubling bifurcation. We also show that the non-autonomous Julian tremor model with a step like increase and decrease in input pressure is able to exhibit oscillations of varying amplitude while keeping a constant frequency spectrum. Both phenomena were observed at Shinmoedake. We also demonstrate that the non-autonomous Julian tremor model with a transient input pressure is able to exhibit long period-like events in addition to harmonic tremor-like events, implying that the same nonlinear mechanism could be responsible for the generation of both type of events.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We present a joint analysis of newly acquired gravity and teleseismic data in the North Tanzanian Divergence, where the lithospheric break-up is at its earliest stage. The impact of a mantle upwelling in more mature branches of the East African Rift has been extensively studied at a lithospheric scale. However, few studies have been completed that relate the deep-seated mantle anomaly detected in broad regional seismic tomography with the surface deformation observed in the thick Archaean Pan-African suture zone located in North Tanzania. Our joint inversion closes the gap between local and regional geophysical studies, providing velocity and density structures from the surface down to ca. 250 km depth with new details. Our results support the idea of a broad mantle upwelling rising up to the lithosphere and creating a thermal modification along its path. However, our study clearly presents an increasing amplitude of the associated anomaly both in velocity and density above 200 km depth, which cannot be solely explained by a temperature rise. We infer from our images the combined impact of melt (2-3%), composition and hydration that accompany the modification of a thick heterogenous cratonic lithosphere are a response to the hot mantle rising. The detailed images we obtained in density and velocity assert that Archaean and Proterozoic units interact with the mantle upwelling to restrict the lithosphere modifications within the Magadi-Natron-Manyara rift arm. The composition and hydration variations associated with those units equilibrate the thermal erosion of the craton root and allow for its stability between 100 and 200 km depth. Above 80 km depth, the crustal part is strongly affected by intruding bodies (melt and gas) which produces large negative anomalies in both velocity and density beneath the main magmatic centers. In addition to the global impact of a superplume, the velocity and density anomaly pattern suggests a 3D distribution of the crust and mantle lithospheric stretching, which is likely to be controlled by inherited fabrics and enhanced by lateral compositional and hydration variations at the Tanzanian craton-orogenic belt boundary.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We report palaeomagnetic and K–Ar geochronologic results of two volcanic sequences from Ethiopia. The Belessa section, dated around 29–30 Ma and spanning ∼1 km in thickness, is related to the Oligocene Afro-Arabian traps, whereas the ∼700-m-thick Debre Sina section was emplaced during the Miocene in two periods around 10–11 and 14–15 Ma. We sampled 67 flows of predominantly basaltic rocks near Belessa and 59 rhyolitic to trachybasaltic flows near Debre Sina. From a geodynamic viewpoint, the magnetostratigraphy of the Belessa sequence confirms that the Ethiopian traps were emplaced at a minimum rate of ∼1 m/kyr, with a possible acceleration of the volume of volcanism over time. To provide insight into the evolution of the geomagnetic field in the Afro-Arabian region over the past 30 Myr, we combined our results with previous studies in the same area. Recentred directional distributions were elongated in the meridian plane, in coherence with field models for a dipole-dominated field. The dispersion 〈span〉S〈/span〉 of the virtual geomagnetic poles, representative of the vigour of palaeosecular variation, was approximately 50% higher during the 10–30 Ma interval than during the past 5 Myr. As the reversal frequency 〈span〉f〈/span〉 was more than two times lower during the Early Oligocene than during the Plio-Pleistocene, it appears that 〈span〉S〈/span〉 and 〈span〉f〈/span〉 are uncorrelated in this near-equatorial region. It remains an open question whether this apparent decoupling is ascribable to a local anomaly, is only sporadic in time, or represents a general feature of the geodynamo.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The non-linear interaction of ocean surface waves produces coherent infrasound noise—microbaroms—between 0.1 and 0.5 Hz. Microbaroms propagate through the atmosphere over thousands of kilometres due to low absorption and efficient ducting between the ground and the stratopause. These signals are globally and permanently detected by the International Monitoring System (IMS) infrasound network, which has been established to monitor compliance with the Comprehensive Nuclear-Test-Ban Treaty. At the International Data Centre (IDC) in Vienna, where IMS data are routinely processed, microbarom detections appear in overlapping frequency bands, and are treated as false alarms. Therefore, understanding the variability in microbarom detections is essential to support the IDC in the reduction of the false alarm rate. In this study, microbarom amplitudes and the direction of arrivals at the German infrasound station IS26 were modelled. For the simulations, the source was described by an operational ocean wave interaction model, and the signal amplitude was modelled using a semi-empirical attenuation relation. This relation strongly depends on middle atmosphere (MA; i.e. 15–90 km altitude) dynamics; however, vertical temperature and wind profiles, provided by numerical weather prediction (NWP) models, exhibit significant biases and differences when compared with high-resolution lidar soundings in altitudes where infrasound signals propagate. To estimate uncertainties in the modelled amplitude, a fully autonomous lidar for MA temperature measurements was installed at IS26. Temperature and wind perturbations, considering observed biases and deviations, were added to the operational high-resolution atmospheric model analysis produced by the European Centre for Medium-Range Weather Forecasts. Such uncertainties in horizontal winds and temperature strongly impact propagation conditions, explaining almost 97 per cent of the actual detections, compared to 77 per cent when using the direct output of the NWP model only. Incorporating realistic wind and temperature uncertainties in NWP models can thus significantly improve the understanding of microbarom detections as well as the detection capability of a single station throughout the year.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The interaction of subducted oceanic lithosphere with the discontinuities of the mantle transition zone (MTZ) provides insight into the composition and temperature of the subducted slab as well as potential melting of the slab or the surrounding mantle and loss of volatiles from the slab. Detailed mapping of the structure of the MTZ will help to better understand how slabs transport material and volatiles into the mantle and how phase transitions affect the slab dynamics. Here we use a dense network of seismic stations in northern Anatolia to image the structure of the MTZ discontinuities in detail using P-wave receiver functions. With a station spacing of about 7 km and a surface footprint of ∼35 km by ∼70 km, analysing receiver functions calculated from teleseismic earthquakes that occurred during an ∼18 month deployment produced clear images of where the mantle transition zone interacts with the Tethys/Cyprus slabs that either lie flat on the 660-km discontinuity or pass into the lower mantle. We observe an undulating 660-km discontinuity depressed by up to 30 km and a slightly depressed (1–2 km) 410-km discontinuity, apparently undisturbed by the slab. The MTZ is thickened to ∼270 km as result of the cool slab in the MTZ influencing the 660-km discontinuity and includes an arrival at ∼520-km depth likely from the top of a flat lying slab or a discontinuity related to a solid-solid phase transition in the olivine component of the mantle. We find evidence for low-velocity zones both above and below the 410-km discontinuity and above the 660-km discontinuity. The low velocity zones around the 410-km discontinuity might be the result of hydration of the MTZ from the slab and upward convection of MTZ material into the upper mantle. The origin of the low velocity zone around the 660-km discontinuity is less clear and could be related to sedimentation of subducted mid-ocean ridge basalts. The small footprint of the seismic array provides accurate information on the structure of the MTZ in an area influenced by subduction and shows small-scale changes in MTZ structure that might be lost in studies covering larger areas with sparser sampling.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Long-period Rayleigh wave horizontal to vertical amplitude (H/V) ratios at a station provide information about local earth structure that is complementary to phase velocity. However, a number of studies have observed that significant scatter appears in these measurements making it difficult to use H/V ratio measurements to resolve earth structure. Some of the scatter in these measurements has been attributed to local geological structure while some has remained unaccounted for. Most Global Seismographic Network (GSN) stations contain two nearby high-quality broad-band seismometers (e.g. in the same vault, but on different piers or in different boreholes). For each broad-band sensor in the IRIS/USGS component of the GSN, we estimate H/V ratios of fundamental mode Rayleigh waves using 〈span〉M〈/span〉 〉 6.5 earthquakes from 2001 to 2018 (around 19 000 measurements). We compute these ratios at a number of discrete periods (25, 50, 75, 100 and 150 s) and find that for well-isolated Rayleigh waves (windows where the correlation coefficients between radial and the phase-shifted vertical components are greater than 0.9) significant scatter in H/V ratios occurs between colocated sensors (greater than 25 per cent at 100 s period). This suggests the scatter in H/V ratio measurements can be at least partially attributed to extremely local phenomena such as sensor emplacement in the vault. We also find that H/V ratios can vary as a function of event backazimuth, indicating that care must be taken when computing average ratios for a station, as a large number of events from a given region could bias H/V ratio measurements at a station.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Full-waveform inversion (FWI) of surface waves can provide a valuable contribution to near-surface investigations, since they are mainly sensitive to the S-wave velocity and hold a high signal-to-noise ratio. We investigate the performance of the individual 2-D elastic FWI of Rayleigh and Love waves as well as the feasibility of a simultaneous joint FWI of both wave types. We apply these three methods to synthetic data as well as to field data. In synthetic reconstruction tests we compare the performance of the individual wave type inversions and explore the benefits of a simultaneous joint inversion. In these tests both individual wave type inversions perform similarly well, given that the initial P-wave velocity model is accurate enough. In the case of an inaccurate initial P-wave velocity model, we observe artifacts in the reconstructed S-wave velocity models of the Rayleigh wave FWI and the joint FWI. If the P-wave velocity is sufficiently known the joint FWI can further improve the model of the S-wave velocity. In the field data application the Love wave FWI is superior to the Rayleigh wave FWI, possibly due to the insufficient initial P-wave velocity model, while the joint FWI further improves the inversion result also in this case.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉As soon as an earthquake starts, the rupture and the propagation of seismic waves redistribute masses within the Earth. This mass redistribution generates in turn a long-range perturbation of the Earth gravitational field, which can be recorded before the arrival of the direct seismic waves. The recent first observations of such early signals motivate the use of the normal mode theory to model the elastogravity perturbations recorded by a ground-coupled seismometer or gravimeter. Complete modeling by normal mode summation is challenging due to the very large difference in amplitude between the prompt elastogravity signals and the direct P-wave signal. We overcome this problem by introducing a two-step simulation approach. The normal mode approach enables a fast computation of elastogravity signals in layered self-gravitating Earth models. The fast and accurate computation of gravity perturbations indicates instrument locations where signal detection may be achieved, and may prove useful in the implementation of a gravity-based earthquake early warning system.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Isostatic gravity anomalies provide a measure of the Earth's gravity field free from the gravitational attractions of the topography and its isostatic compensation, most commonly represented by a variation in the depth of a compensating density contrast, for example the Moho. They are used by both geodesists and geophysicists alike, though often for different purposes. Unfortunately though, the effect of subsurface loading on the lithosphere renders transfer function (admittance) methods unusable when surface and subsurface loads coexist. Where they exist, subsurface loads are often expressed in the Bouguer anomaly but not in the topography, and it is shown here that this phase disconnect cannot be faithfully represented by either real- or complex-valued analytic admittance functions. Additionally, many studies that employ the isostatic anomaly ignore the effects of the flexural rigidity of the lithosphere, most often represented as an effective elastic thickness (〈span〉Te〈/span〉), and assume only Airy isostasy, i.e. surface loading of a plate with zero elastic thickness. The consequences of such an omission are studied here, finding that failure to account for flexural rigidity and subsurface loading can result in (1) over- or underestimates of both inverted Moho depths and dynamic topography amplitude, and (2) underestimates of the size of topographic load that can be supported by the plate without flexure. An example of the latter is shown over Europe. Finally, it is demonstrated how low values of the isostatic anomaly variance can actually be biased by these anomalies having low power at the long wavelengths while still possessing high power at middle to short wavelengths, compared to the corresponding Bouguer anomaly power spectrum. This will influence the choice of best-fitting isostatic model if the model is chosen by minimisation of the isostatic anomaly standard deviation.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Anisotropic parameters, when considered in a microseismic processing, are typically inverted using perforation shot data and/or simultaneously inverted with microseismic event locations. Inverting using microseismic data alone usually leads to an under-constrained inverse problem that is highly dependent on prior/initial information. We carefully processed the waveforms from perforation shots, and picked P-, SH-, and SV-wave arrival times to mitigate this issue. Because the perforation shot locations are known, the inversion is better-constrained by reducing the number of model parameters while increasing the number of observations. We applied both Maximum A Posteriori (MAP) estimation and Randomized Maximum Likelihood (RML) methods for anisotropic parameter estimation, uncertainty quantification, and trade-off analysis. Results verified the stability of the inversion and revealed the uncertainty and trade-offs among model parameter. In addition, attenuation is generally not considered in microseismic modeling and processing. Our study found that hydraulic stimulation may lead strong increases in seismic attenuation to reservoirs. The attenuation can dramatically change waveform characteristics and cause velocity dispersion. Thus, sonic logs, which are acquired at frequencies much higher than seismic data frequency should not be used directly for data processing in hydraulic stimulated zones.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉On 26 May 2006 at 23:54 UTC, a moderate shallow crustal earthquake with a moment magnitude of 6.3 occurred in the southern part of Yogyakarta in Java, Indonesia. The earthquake caused severe damages in the area in addition to over 5700 casualties. The cause of this earthquake was initially believed to have been a rupture on the northeast-southwest trending Opak Fault; however, the role of this fault in the earthquake continues to be debated. Therefore, this study presents a subsurface model constructed to characterize the fault geometry associated with the earthquake. We used previously reported aftershock data to image subsurface velocity variations through seismic tomographic inversion of primary waves, shear waves, and their velocity ratio (V〈sub〉p〈/sub〉/V〈sub〉s〈/sub〉). Using data from 10 stations around the hypothetical fault, 588 aftershock events were mostly located 10–15 km east of the Opak River Fault with a maximum depth of approximately 20 km. The seismic tomographic inversion results indicated that severe damage during the earthquake occurred in areas with larger V〈sub〉p〈/sub〉/V〈sub〉s〈/sub〉 ratios associated with unconsolidated sediments, in accordance with previous findings. Furthermore, the configuration of an unnamed fault that was activated during the earthquake is delineated by a velocity anomaly with a depth of up to 5–7 km. This structure is interpreted as a strike slip fault with a reverse component dipping to the east, striking northeast-southwest.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Major faults release tectonic stress not only during large earthquakes, but also through slow aseismic slip associated with microseismic activity. Micro earthquakes could be interpreted as the failure of fault asperities. If the slip accumulated during microearthquakes could be estimated with seismological data, the aseismic component of fault slip is more difficult to evaluate. In this study, we are interested in quantifying the aseismic fault slip from the seismic moment released by the simultaneous on fault seismic activity. For that, we study systematically the relationship between total slip and microseismic moment released on a fault modeled as a heterogeneous planar Dieterich-Ruina frictional interface producing earthquakes and aseismic slip. It is shown that for a long enough earthquake sequence, typically involving several ocurrences of the largest asperity failure, the seismic moment released is directly proportional to the total (seismic plus aseismic) slip on the fault accumulated during the sequence. In particular, the total seismic moment divided by the fault area represents a fraction of the total slip that is of the order of the apparent density of asperities. The apparent density of asperities is defined as the fraction of the fault area that experienced a dynamic rupture during the sequence. This result provides a strong mechanical support to extend the use of microseismic activity as creepmeters along major active faults, beyond the common analysis of repeating earthquakes.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We introduce in the on-site earthquake early warning (EEW) a partially non-ergodic perspective from the site effects point of view. We consider the on-site EEW approach where the peak ground velocity (PGV) for S-waves is predicted from an early estimate, over the P-waves, of either the peak-displacement (PD) or cumulative squared velocity (IV2). The empirical PD-PGV and IV2-PGV relationships are developed by applying a mixed-effect regression where the site-specific modifications of ground shaking are treated as random effects. We considered a large data set composed of almost 31000 selected recordings in central Italy, a region struck by four earthquakes with magnitude between 6 and 6.5 since the 2009 L’Aquila earthquake. We split the data set into three subsets used for calibrating and validating the on-site EEW models, and for exemplifying their application to stations installed after the calibration phase. We show that the partially non-ergodic models improve the accuracy of the PGV predictions with respect to ergodic models derived for other regions of the world. Moreover, considering PD and accounting for site effects, we reduce the (apparent) aleatory variability of the logarithm of PGV from 0.31-0.36, typical values for ergodic on-site EEW models, to about 0.25. Interestingly, a lower variability of 0.15 is obtained by considering IV2 as proxy, which suggests further consideration of this parameter for the design of on-site EEW systems. Since being site-specific is an inherent characteristic of on-site EEW applications, the improved accuracy and precision of the PGV predicted for a target protection translate in a better customization of the alert protocols for automatic actions.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We calculate the change in effusion rate of lava from a volcanic fissure due to pressure changes in the volcanic conduit. The conduit is modelled as a cylinder with elliptical cross section, embedded in an elastic medium. The elliptical shape can represent a wide range of cross sections, according to the value of eccentricity, from almost circular vents to very long and narrow fissures. A 2-D problem is considered assuming invariance of pressure changes and conduit geometry with depth. The problem is solved analytically and expressions for the displacement and the stress fields in the elastic medium are provided. The displacement of the conduit wall is proportional to the ratio between the pressure change and the rigidity of surrounding rocks. The flow rate is a nonlinear function of the pressure change and increases with increasing pressure, due to the elastic deformation of the conduit wall. We consider flow rate oscillations with periods ranging from several minutes to hours, as are often observed during effusive eruptions. Assuming pressure oscillations with these periods, flow rate oscillations resulting from the elastic deformation of the conduit are calculated. The greatest oscillations in flow rate are obtained for very large values of the conduit eccentricity, corresponding to long and narrow volcanic fissures. For example, if a fissure is 100 m long and 2 m large, a pressure oscillation with an amplitude of 1 MPa yields a maximum displacement of the conduit wall equal to about 6 cm and an amplitude of flow rate oscillations of about 20%.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We propose a numerical method to compute the inertial modes of a container with near-spherical geometry based on the fully spectral discretisation of the angular and radial directions using spherical harmonics and Gegenbauer polynomial expansion respectively. This allows to solve simultaneously the Poincaré equation and the no penetration condition as an algebraic polynomial eigenvalue problem. The inertial modes of an exact oblate spheroid are recovered to machine precision using an appropriate set of spheroidal coordinates. We show how other boundaries that deviate slightly from a sphere can be accommodated for with the technique of 〈span〉equivalent spherical boundary〈/span〉 and we demonstrate the convergence properties of this approach for the triaxial ellipsoid.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We outline a method using Gradient Flow Independent Component Analysis (ICA) to separate signals comprising the coda in a topographically complex setting. We also identify the sources of scattered signals by tracking signal backazimuths over time. The Gradient Flow ICA method is shown to effectively separate signals in the acoustic coda. The method correctly identifies the backazimuth of the first arrival from two 800 kg TNT equivalent explosions as well as subsequent signals scattered by the surrounding topography. Circular statistics is used to analyze the variance, mean, and uniformity of calculated backazimuths. These results have strong implications for understanding the acoustic wavefield by identifying scatterers and inverting for atmospheric conditions.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉A detailed and precise knowledge of ocean bottom topography is essential in many geoscientific and oceanographic applications. Shipborne echo sounding provides the only direct bathymetric method. However, even after decades of applying this technique only a fraction of the global ocean could be covered. Alternatively, gravity data inversion is a feasible method to infer ocean bottom topography since the gravity field correlates with topography at short to medium wavelengths. Gravity field observables are globally provided by dedicated satellite missions like GOCE, GRACE and GRACE-FO and, over the oceans, by satellite altimetry. Regional and local measurements are realized by means of ground-based, shipborne and airborne gravimetry. For the first time in Europe a jet aircraft was used for airborne gravimetry. During the GEOHALO flight campaign over Italy and the Tyrrhenian, Ionian and Adriatic Seas in June 2012, the German research aircraft HALO carried an entire suite of geodetic-geophysical instrumentation, including gravity meters. The careful processing of the gravity data acquired by the CHEKAN-AM instrument of the German Research Centre of Geosciences allowed to achieve an accuracy at the mGal level. Subsequently, the Parker-Oldenburg inversion was applied to predict ocean bottom topography along the GEOHALO profiles flown over the ocean. To constrain the parameter space in the inversion we used the Crust1.0 model. Finally, the obtained results were compared to the General Bathymetric Chart of the Oceans (GEBCO) in order to estimate the performance of the method. Our study demonstrates that airborne gravimetry aboard a jet aircraft is capable to provide valuable data for regional geoscientific studies.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉A new method is proposed to estimate separately 3D heterogeneities of intrinsic and scattering attenuation parameters of seismic shear waves. The method, which assumes isotropic scattering, consists of two steps. First, the source, site, and mean path (both intrinsic and scattering attenuation) terms are estimated by an envelope-fitting method for each earthquake event; second, these attenuation parameters are mapped to 3D space in a tomographic manner by using sensitivity kernels of the attenuation parameters. Monte Carlo simulation is adopted to take account of the effect of the depth-dependent velocity structure in the envelope calculation. The proposed method was then applied in south-western Japan. The result showed that both intrinsic and scattering attenuation were stronger in the Kyushu region, where tectonic activity was higher, than in other studied regions. This association of attenuation parameters with tectonic activity is consistent with previous findings. The proposed method thus provides an alternative to multiple lapse time window analysis, a widely used method for separately estimating intrinsic and scattering attenuation factors.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We discuss a regularization approach to improving numerical stability of the frequency domain Maxwell equations of electromagnetic geophysics. To enforce divergence-free conditions we add a scaled grad-div operator to the curl-curl equation for electric fields. This deflates the null space of the curl-curl operator and significantly reduces the condition number of the linear system derived for the finite difference (FD) approximation, resulting in faster and more stable iterative solution. We explicitly discuss extensions needed to solve the adjoint problems and consider two apparently different approaches to the sensitivity calculation, which we show are ultimately equivalent. To complete assessment of the new approach in practice we have implemented the modified solver in the ModEM inversion code, and compared it to more standard methods (based on solving with a divergence correction) on inversion of a real data set. This, and simpler tests of the forward solver in isolation, demonstrate that, with appropriate scaling of the added grad-div term, the regularization approach can dramatically improve the efficiency of Krylov subspace methods. Run times for the inversion test are reduced by almost a factor of three, making 3D FD inversion significantly more practical.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Anisotropic viscosity is often discussed in the context of flow-induced crystallographic alignment, also known as lattice-preferred orientation, believed to exist throughout the mantle. However, other geophysical systems can produce bulk anisotropic rheology. One example is so-called shape-preferred orientation, such as that imparted by repeated, non-randomly oriented igneous intrusions. Dikes and sills tend to intrude perpendicular to least compressive tectonic stress, and therefore are typically aligned with other coeval intrusions. Intrusions provide distinct planes of weakness, producing, in a bulk sense, an anisotropic fabric to the rocks in which they are intruded. We parameterize this style of anisotropic viscosity with an effective-medium approximation, based on analytically derived characteristics of layered media. We find that the magnitude of anisotropy is highly sensitive to the presence of weak layers, even if they are relatively thin and/or sparsely distributed. We test the effects of our constitutive law on Rayleigh-Taylor instability growth, representing lithospheric instabilities. Our results indicate dramatically decreased stability of lithosphere intruded by dikes and/or sills.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We use one month of continuous seismic waveforms from a very dense seismic network to image with unprecedented resolution the shallow damage structure of the San Jacinto fault zone across the Clark fault strand. After calculating noise correlations, high apparent-velocity arrivals coming from below the array are removed using a frequency-wavenumber filter. This is followed by a double-beamforming analysis on multiple pairs of subarrays to extract phase and group velocity information across the study area. The phase and group velocity dispersion curves are regionalized into phase and group velocity maps at different frequencies, and these maps are inverted using a neighbourhood algorithm to build a 3D shear wave velocity model around the Clark fault down to ∼500 m depth. The model reveals strong lateral variations across the fault strike with pronounced low velocity zones corresponding to a local sedimentary basin and a fault zone trapping structure. The results complement previous earthquake- and seismic noise-based imaging of the fault zone at greater depth and clarify properties of structural features near the surface.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉This paper studies the effect of omission error on height system unification by analyzing the contributions of the low and high spherical harmonic degree components of a local vertical datum (LVD). Previous studies have combined the satellite-only GOCE (Gravity field and steady-state Ocean Circulation Explorer) global geopotential models (GGMs) with a high-resolution GGM (and in some cases topography-implied gravity field signals) in order to account for (or minimize) the omission error of GOCE. This study attempts to answer the question of whether the omission error should be accounted for and to what extent by decomposing the Brazilian Vertical Datum referenced to a tide gauge station at Imbituba (BVD-I) into low and high degree components using the Gaussian averaging function. A combined EIGEN-GOCE GGM is used to compute the high degree datum offset. Overall, the results confirm that (1) the low degree component accounts for the majority of the discrepancy between the BVD-I and the global datum but is insufficient to accurately determine the BVD-I offset parameter, (2) the high degree component contributes little to the overall discrepancy between the datums, but significantly improves the errors in the low degree component, and (3) the signals beyond the resolution of the high-resolution GGM do not significantly impact the datum offset parameters, so that it is not necessary to account for the omission error of the high-resolution GGM using topography-implied gravity field signals. The approach used in this work adds little computational cost, especially over large spatial extents and dense GPS/leveling data, compared to the computation of topography-generated gravity field quantities and has the added advantage of being able to determine the low degree systematic errors in the GPS/leveling data (LVD) using the unbiased satellite-only GOCE models as reference.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Earth rotation studies require accurate knowledge of the global oceanic velocity and mass fields, for proper accounting of ocean angular momentum (OAM) effects on the planetary budget. We analyze a new OAM series (1992–2015) based on the solution of a global general circulation model constrained to most existing ocean data. The impact of the data-constrained optimization on OAM is substantial, and particularly essential for calculating effects of global mean ocean mass changes, which can be important for determining annual cycles and long term trends in OAM. The contributions of sea ice to OAM variations, also estimated, are found to be negligible. Uncertainties in OAM series are assessed by comparison with other available estimates. Results indicate low signal-to-noise ratios for all the analyzed OAM series. Comparisons with geodetic, atmospheric and hydrologic data, in the context of the planetary angular momentum budget, point to the continued need for improvements in some or all of the series. Possible paths are offered for producing better OAM estimates in the future.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We present 96 new seafloor heat flow determinations, made with a 3.5-m violin-bow probe and collocated with seismic reflection profiles, from the northern Hikurangi margin on the east coast of New Zealand's North Island. Here the Hikurangi Plateau on the Pacific Plate is subducting under the Australian Plate. The background heat flow is 58 ± 8 mW m〈sup〉−2〈/sup〉, consistent with and within the variability of globally observed heat flow (56 ± 15 mW m〈sup〉−2〈/sup〉) for oceanic crust 90 to 120 Ma, the age of the Hikurangi Plateau. Seaward of the deformation front, we find evidence for advective fluid flow associated with basement relief. Landward of the deformation front, we use a two-dimensional steady-state finite element model to quantify the thermal regime. Despite corrections for the effects of bottom water temperature change, bathymetry, and sedimentation there is considerable scatter in the heat flow data including a local and sharp increase in heat flow of up to 35 mW m〈sup〉−2〈/sup〉 observed over the outermost wedge. Variability in the heat flow data is likely due to complex, unmodeled three-dimensional fluid flow. We augment our heat flow measurements with estimates from a bottom-simulating reflection (BSR) and continental bottom hole temperatures and conclude that the effective coefficient of friction,  μ*, is approximately 0.06 in the region of observed slow slip events and increases to 0.18 approximately 50 km landward of the deformation front. This transition in μ* may be marking the downdip edge of overpressures along the subduction thrust, suggesting that slow slip is enabled by overpressure.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We analyze noise characteristics of spatial and temporal correlation of 260 continuous GPS sites from Crustal Movement Observation Network of China (CMONOC). These datasets were mainly collected between 2010 and 2016, with an average of 6 years of position time series. In the functional analysis, a clear regional dependence of seasonal movements has been observed and other significant periodic signals are detected nearby the GPS draconitic period and its harmonics. The distribution of these periodic signals shows a spatial correlation, along with non-negligible local inconsistencies. In the stochastic analysis, impacts of the periodicities on the noise assessment have been investigated and Maximum likelihood estimation is used to study noise properties for deseasonalized residual time series having the seasonal signals removed and filtered residual time series having other periodic signals removed further. We demonstrate that for both solutions, the flicker noise is thought to be the dominant time-correlated noise and velocity uncertainties may be underestimated 8–10 times if assuming a pure white noise. Ignoring the periodicities could bias the estimation of noise amplitude, spectral index and velocity uncertainty. After removing the periodic signals, the median flicker noise magnitude shows an average of ∼10 per cent reductions and the noise process shifts closer toward white noise. A correlation between the index variations and RMS variations has been observed, indicating that the index varies more significantly for sites with more periodic signals removed. Besides, the spectral index in the vertical component has a better spatial correlation than that in the horizontal and the spatial distribution of the index of deseasonalized solutions seems to correlate well with the amplitudes of seasonal signals, probably implying common sources of spatial variations of these characteristics. Furthermore, analyses of inter-site correlations indicate that the correlation induced from the periodic signals displays a similar pattern to the deseasonalized solutions, confirming that the period pattern is spatially correlated and can induce time-series correlations. Finally, the stochastic processes of the common mode noise are predominantly featured by spatially correlated flicker and white noise over a wide range, consistent with previous results.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Waveform backprojection is a key technique of earthquake-source imaging, which has been widely used for extracting information of earthquake source evolution that cannot be obtained by kinematic source inversion. The technique enjoys considerable popularity, owing to the simplicity of its implementation and the robustness of its processing, but the physical meaning of backprojection images has remained elusive. In this study, we reviewed the mathematical representation of backprojection (BP) and hybrid backprojection (HBP) methods, following the pioneering work of Fukahata et al. (Geophys. J. Int. (2014) 196, 552559), to clarify the physical implications of BP images. We found that signal intensity in BP and HBP images is scaled with the amplitude of the Green’s function that corresponds to a unit-step slip, which results in the signal intensity being depth dependent. We propose variants of BP and HBP, which we call kinematic BP and HBP, respectively, to relate the BP signal intensity to slip motion of an earthquake by modifying the normalizing factors used in the original BP and HBP methods. The original BP and HBP images remain useful for assessing the spatiotemporal strength of the wave radiation, which scales with the amplitude of the Green’s function, whereas the kinematic BP and HBP methods are suitable for imaging the slip motion that is responsible for the high-frequency radiation produced during the source-rupture process.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉In the field of controlled source seismology, the acoustic 3-D Full Waveform Inversion (FWI) technique has become a common tool for imaging geologically complex structures in land, as well as in marine settings. However, the Earth behaves elastically and therefore, excluding the elastic effect could have a significant impact on the inversion results, especially for interfaces with large S-wave velocity contrasts. To examine the contribution of the elastic approach, we compare acoustic and elastic 3-D FWI applied to a 3-D seismic dataset from the East Pacific Rise (EPR) 9°50′ N, where the subsurface is represented by igneous basaltic rocks. To establish an efficient inversion strategy, we first conducted a number of tests, which suggest a simultaneous, multi-parameter inversion as the most efficient approach when inverting signals with frequencies below 7 Hz. The reduction in the total misfit for the elastic case is 10–15 per cent lower than that for the acoustic one, suggesting that the elastic approach explains the observed data better than the acoustic approach. Furthermore, the compressional velocity images of the upper-oceanic crust obtained using the two approaches differ significantly, not only in velocity magnitude but also structurally. We argue that the results obtained from the acoustic modeling are geologically less plausible and suggest the elastic model as a more reliable representation of the upper-oceanic crust.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We investigated an Mw∼6.2 earthquake doublet on the border of the USA and Canada using ALOS2 Light-of-Sight displacements and GPS measurements. We selected three L-band ALOS-2 interfergorams with temporal baselines of one year to extract coseismic deformation maps, in which master and slave images were both acquired in July. A subpixel-based alignment and range spectral splitting techniques under the GAMMA InSAR software framework were applied to improve the interferometric coherence and reduce the effects of phase anomalies in two of the three interferometric pairs due to either ionospheric delay or a potential focusing issues in the generation of the ALOS2 SLC data. The updated interferograms convincingly reveal deformation fringe patterns produced by the two earthquakes. We conducted a nonlinear geophysical inversion to estimate the geometric parameters of the earthquakes with the InSAR and GPS measurements. The best fitting model shows that a thrust faulting on a reverse fault and left-lateral strike-slip faulting on a nearly vertical fault with the centroid depths of 6.7 and 8.5 km, respectively, are most likely responsible for the earthquake doublet. The eastern Denali fault (EDF) and Duke River fault are major active faults in the region and the earthquake doublet could be due to reactivation of the part of the two faults system.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉We present the theory for and applications of Hamiltonian Monte Carlo (HMC) solutions of linear and nonlinear tomographic problems. HMC rests on the construction of an artificial Hamiltonian system where a model is treated as a high-dimensional particle moving along a trajectory in an extended model space. Using derivatives of the forward equations, HMC is able to make long-distance moves from the current towards a new independent model, thereby promoting model independence, while maintaining high acceptance rates. Following a brief introduction to HMC using common geophysical terminology, we study linear (tomographic) problems. Though these may not be the main target of Monte Carlo methods, they provide valuable insight into the geometry and the tuning of HMC, including the design of suitable mass matrices, and the length of Hamiltonian trajectories. This is complemented by a self-contained proof of the HMC algorithm in the Appendix. A series of tomographic/imaging examples is intended to illustrate (i) different variants of HMC, such as constrained and tempered sampling, (ii) the independence of samples produced by the HMC algorithm, and (iii) the effects of tuning on the number of samples required to achieve practically useful convergence. Most importantly, we demonstrate the combination of HMC with adjoint techniques. This allows us to solve a fully nonlinear, probabilistic traveltime tomography with several thousand unknowns on a standard laptop computer, without any need for supercomputing resources.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉 〈strong〉Moment Tensors of Hydraulically-Induced AEs:〈/strong〉 Hydraulic fracturing is an important technique in the development of enhanced geothermal systems and unconventional resources. Although the fracture modes induced by hydraulic fracturing influence the recovery efficiency of the resources, the current understanding of this relationship is insufficient. In the present study, we considered the acoustic emissions (AEs) induced during hydraulic fracturing under uniaxial loading conditions in the laboratory, and applied a moment tensor analysis by carefully correcting the coupling condition and directivity of AE transducers. Experiments were conducted for two types of Kurokami-jima granite samples: those with a rift plane perpendicular (Type H) or parallel (Type V) to the expected direction of fracture propagation (i.e. along the loading axis). In the experiments, both sample types experienced a significant number of shear, tensile, and compressive events. The dominant fracture mode for Type H samples is found to be tensile events in which the fracture plane is parallel to the loading axis, whereas for Type V samples, shear events are dominant. This difference suggests that the dominant fracture modes induced by hydraulic fracturing are highly dependent on the relationship between the direction of fracture propagation and orientation of pre-existing weak planes.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉The size of an earthquake can be defined either from the seismic moment (M〈sub〉0〈/sub〉) or in terms of radiated seismic energy (E〈sub〉r〈/sub〉). These two parameters look at the source complexity from different perspectives: M〈sub〉0〈/sub〉 is a static measure of the earthquake size, whereas E〈sub〉r〈/sub〉 is related to the rupture kinematics and dynamics. For practical applications and for dissemination purposes, the logarithms of M〈sub〉0〈/sub〉 and E〈sub〉r〈/sub〉 are used to define the moment magnitude M〈sub〉w〈/sub〉 and the energy magnitude M〈sub〉E〈/sub〉, respectively. The introduction of M〈sub〉w〈/sub〉 and M〈sub〉E〈/sub〉 partially obscure the complementarity of M〈sub〉0〈/sub〉 and E〈sub〉r〈/sub〉. The reason is due to the assumptions needed to define any magnitude scale. For example, in defining M〈sub〉w〈/sub〉, the apparent stress (i.e. the ratio between M〈sub〉0〈/sub〉 and E〈sub〉r〈/sub〉 multiplied by the rigidity) was assumed to be constant, and under this condition, M〈sub〉w〈/sub〉 and M〈sub〉E〈/sub〉 values would only differ by an off-set which, in turn, depends on the average apparent stress of the analysed dataset. In any case, when the apparent stress is variable and, for example, scales with M〈sub〉0〈/sub〉, the value of M〈sub〉E〈/sub〉 derived from M〈sub〉w〈/sub〉 cannot be used to infer E〈sub〉r〈/sub〉.In this study, we investigate the similarities and differences between M〈sub〉w〈/sub〉 and M〈sub〉E〈/sub〉 in connection with the scaling of the source parameters using a dataset of around 4700 earthquakes recorded at both global and regional scales and belonging to four datasets. These cover different geographical areas and extensions and are composed by either natural or induced earthquakes in the magnitude range 1.5 ≤ M〈sub〉w〈/sub〉 ≤ 9.0. Our results show that M〈sub〉E〈/sub〉 is better than M〈sub〉w〈/sub〉 in capturing the high-frequency ground shaking variability whenever the stress drop differs from the reference value adopted to define M〈sub〉w〈/sub〉. We show that M〈sub〉E〈/sub〉 accounts for variations in the rupture processes, introducing systematic event-dependent deviations from the mean regional peak ground motion velocity scaling. Therefore, M〈sub〉E〈/sub〉 might be a valid alternative to M〈sub〉w〈/sub〉 for deriving ground motion prediction equations for seismic hazard studies in areas where strong systematic stress drop scaling with M〈sub〉0〈/sub〉 are found, such as observed for induced earthquakes in geothermal regions. Furthermore, we analyse the different datasets in terms of their cumulative frequency-magnitude (CFM) distribution, considering both M〈sub〉E〈/sub〉 and M〈sub〉w〈/sub〉. We show that the b values from M〈sub〉w〈/sub〉 (〈span〉b〈/span〉〈sub〉Mw〈/sub〉) and M〈sub〉E〈/sub〉 (〈span〉b〈/span〉〈sub〉ME〈/sub〉) can be significantly different when the stress drop shows a systematic scaling relationship with M〈sub〉0〈/sub〉. We found that 〈span〉b〈/span〉〈sub〉ME〈/sub〉 is nearly constant for all datasets, while 〈span〉b〈/span〉〈sub〉Mw〈/sub〉 shows an inverse linear scaling with apparent stress.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Changes in the Earth's magnetic field can deeply modify the polar caps and auroral zones, which are the regions of most frequent precipitation of energetic particles. The present field is characterized by a dominant dipole plus weaker multipolar components. The field varies greatly in time, with the most drastic changes being polarity reversals that take place on average every ∼200,000 years. During a polarity transition the field magnitude may diminish to about 10 per cent of its value prior to the reversal due to a decreasing dipolar component and by becoming mostly multipolar in nature. Polar caps depend on the geomagnetic field configuration so changes in their morphology are expected as a consequence of the variation and reversal of this field. We model polar caps’ location by considering a superposition of the internal geomagnetic field and a uniform external field and then following the open field lines to Earth's surface. Polar caps’ location and shape for different magnetic field reversal scenarios are analyzed in the present work. Two polar caps near the present dipole axis intersection with Earth's surface prevail for a dipole decrease to a certain extent, below which the southern hemisphere polar cap moves to mid-latitudes. An axial dipole collapse gives a pair of polar caps both at mid-latitudes of the southern hemisphere, while in a dipole rotation scenario the polar caps reside at the equator. If reversals occur due to an energy cascade from the dipole to higher degrees, more than two polar caps may appear. In our energy cascade scenario four polar caps at various latitudes of both hemispheres prevail. These results indicate that during reversals auroral zones may reach mid and low-latitude regions, and the atmosphere may become more vulnerable to the direct effect of energetic particle precipitation. This vulnerability is particularly striking at the southern hemisphere where reversed flux patches appear on the core-mantle boundary and weak intensity characterizes the present field at Earth's surface.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉It is strongly debated whether the interface between the lithosphere and underlying asthenosphere is a temperature-dependent rheological transition, as expected in a thermal convection system, or additionally affected by the presence of melts and/or fluids. Previous surface wave studies of Pacific oceanic lithosphere have found that shear velocity and azimuthal anisotropy vary with seafloor crustal age as expected for a thermal control; however radial anisotropy does not. Various thermo-mechanical models have been proposed to explain this disparate behaviour. Nonetheless, it is unclear how robust the surface wave constraints are, and this is what we test in this study. We apply a Bayesian model space search approach to three published Pacific surface-wave dispersion datasets, two phase-velocity and one combined phase- and group-velocity set, and determine various proxies for the depth of the lithosphere-asthenosphere boundary (LAB) and their uncertainties based on the velocity and radial anisotropy model distributions obtained. In their overall character and pattern with age, the velocity models from different datasets are consistent with each other, although they differ in their values of LAB depths. Uncertainties are substantial (as much as 20 km on LAB depths) and the addition of group-velocity data does not reduce them. Radial anisotropy structures differ even in pattern and display no obvious age dependence. However, given the uncertainties, we cannot exclude that radial anisotropy, azimuthal anisotropy, and velocity models actually reflect compatible, age-dependent, LAB depth estimates. The velocity LAB trends are most like those expected for half-space cooling, because velocity differences persist at old ages, below the depth of common plate cooling models. Any direct signature of sub-ridge melt would be too small-scale to be resolved by these data. However, the velocity-increasing effects of dehydration and depletion due to melting below the ridge could explain why LAB proxy depths tend to a minimum of ∼60 km below young ocean floor.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Seismic radiation from indigenous sources can be represented by the excess of model stress over actual stress, a second-order tensor field that Backus & Mulcahy (1976) named the 〈span〉stress glut〈/span〉. We prove a new representation theorem that exactly and uniquely decomposes any stress-glut (or strain-glut) density into a set of orthogonal tensor fields of increasing degree $\alpha$, up to six in number, ordered by their first nonzero polynomial moments. The zeroth-degree field ($\alpha \ = \ 0$) is the projection of the stress-glut density onto its zeroth polynomial moment, which defines the seismic moment tensor, Aki seismic moment ${M_0}$, and centroid-moment tensor (CMT) point source. The higher-degree fields describe mechanism complexity—source variability that arises from the spacetime variations in the orientation of the stress glut, which is the main subject of this theoretical study. The representation theorem generalizes the point-source approximation to a sum of multipoles that features the CMT monopole as its leading term. The first-degree field contributes a dipole tensor with an mechanism orthogonal to the CMT, the second-degree field contributes a quadrupole tensor, and so on, up to six orthogonal fields in all. We define the total scalar moment ${M_{\rm{T}}}$ to be the integral of the scalar moment density, and we use the representation theorem to partition this total moment into a sum of fractional moments for each degree $\alpha $. If the faulting is simple enough, ${M_{\rm{T}}} \approx {M_0}$, and the higher-degree terms will be small. When the faulting is more complex, however, ${M_{\rm{T}}} 〉 {M_0}$; the higher-degree fields will contribute more to the radiation, and this contribution will increase with frequency. Application to simple planar faulting shows that out-of-plane variations in slip-vector orientation reduce ${M_0}/{M_{\rm{T}}}$ more than in-plane variations of similar magnitude, consistent with Kagan's (1988) early study of mechanism complexity. We decompose stress-glut realizations from the Graves & Pitarka (2016) rupture simulator; typical values of ${M_0}/{M_{\rm{T}}}$ are 0.82–0.92, consistent with analytical results. The higher-degree fields of the Graves-Pitarka sources typically radiate up to $\alpha \ = \ 4$; only the isotropic term is zero. We compute synthetic seismograms to illustrate the radiation patterns of the higher-degree fields and their frequency dependence. Decomposition of the Duputel & Rivera (2017) source model for the 2016 Kaikoura earthquake (M〈sub〉w〈/sub〉 7.8) indicates that the radiation from the higher-degree fields was large enough (${M_0}/{M_{\rm{T}}}$ = 0.82) that it may be possible to invert global datasets for low-degree multipoles. Other applications of the stress-glut representation theorem are discussed.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We present a joint analysis of newly acquired gravity and teleseismic data in the North Tanzanian Divergence, where the lithospheric break-up is at its earliest stage. The impact of a mantle upwelling in more mature branches of the East African Rift has been extensively studied at a lithospheric scale. However, few studies have been completed that relate the deep-seated mantle anomaly detected in broad regional seismic tomography with the surface deformation observed in the thick Archaean Pan-African suture zone located in North Tanzania. Our joint inversion closes the gap between local and regional geophysical studies, providing velocity and density structures from the surface down to 〈span〉ca〈/span〉. 250  km depth with new details. Our results support the idea of a broad mantle upwelling rising up to the lithosphere and creating a thermal modification along its path. However, our study clearly presents an increasing amplitude of the associated anomaly both in velocity and density above 200 km depth, which cannot be solely explained by a temperature rise. We infer from our images the combined impact of melt (2–3 per cent), composition and hydration that accompany the modification of a thick heterogenous cratonic lithosphere are a response to the hot mantle rising. The detailed images we obtained in density and velocity assert that Archaean and Proterozoic units interact with the mantle upwelling to restrict the lithosphere modifications within the Magadi–Natron–Manyara rift arm. The composition and hydration variations associated with those units equilibrate the thermal erosion of the craton root and allow for its stability between 100 and 200  km depth. Above 80 km depth, the crustal part is strongly affected by intruding bodies (melt and gas) which produces large negative anomalies in both velocity and density beneath the main magmatic centres. In addition to the global impact of a superplume, the velocity and density anomaly pattern suggests a 3-D distribution of the crust and mantle lithospheric stretching, which is likely to be controlled by inherited fabrics and enhanced by lateral compositional and hydration variations at the Tanzanian craton-orogenic belt boundary.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉The nonlinear interaction of ocean surface waves produces coherent infrasound noise—microbaroms—between 0.1 and 0.5 Hz. Microbaroms propagate through the atmosphere over thousands of kilometres due to low absorption and efficient ducting between the ground and the stratopause. These signals are globally and permanently detected by the International Monitoring System (IMS) infrasound network, which has been established to monitor compliance with the Comprehensive Nuclear-Test-Ban Treaty. At the International Data Centre (IDC) in Vienna, where IMS data are routinely processed, microbarom detections appear in overlapping frequency bands, and are treated as false alarms. Therefore, understanding the variability in microbarom detections is essential to support the IDC in the reduction of the false alarm rate. In this study, microbarom amplitudes and the direction of arrivals at the German infrasound station IS26 were modelled. For the simulations, the source was described by an operational ocean wave interaction model, and the signal amplitude was modelled using a semi-empirical attenuation relation. This relation strongly depends on middle atmosphere (MA; i.e. 15–90 km altitude) dynamics; however, vertical temperature and wind profiles, provided by numerical weather prediction (NWP) models, exhibit significant biases and differences when compared with high-resolution light detection and ranging instrument (lidar) soundings in altitudes where infrasound signals propagate. To estimate uncertainties in the modelled amplitude, a fully autonomous lidar for MA temperature measurements was installed at IS26. Temperature and wind perturbations, considering observed biases and deviations, were added to the operational high-resolution atmospheric model analysis produced by the European Centre for Medium-Range Weather Forecasts. Such uncertainties in horizontal winds and temperature strongly impact propagation conditions, explaining almost 97 per cent of the actual detections, compared to 77 per cent when using the direct output of the NWP model only. Incorporating realistic wind and temperature uncertainties in NWP models can thus significantly improve the understanding of microbarom detections as well as the detection capability of a single station throughout the year.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉In the field of controlled source seismology, the acoustic 3-D Full Waveform Inversion (FWI) technique has become a common tool for imaging geologically complex structures in land, as well as in marine settings. However, the Earth behaves elastically and, therefore, excluding the elastic effect could have a significant impact on the inversion results, especially for interfaces with large 〈span〉S〈/span〉-wave velocity contrasts. To examine the contribution of the elastic approach, we compare acoustic and elastic 3-D FWI applied to a 3-D seismic data set from the East Pacific Rise (EPR) 9°50′ N, where the subsurface is represented by igneous basaltic rocks. To establish an efficient inversion strategy, we first conducted a number of tests, which suggest a simultaneous, multiparameter inversion as the most efficient approach when inverting signals with frequencies below 7 Hz. The reduction in the total misfit for the elastic case is 10-15 per cent lower than that for the acoustic one, suggesting that the elastic approach explains the observed data better than the acoustic approach. Furthermore, the compressional velocity images of the upper-oceanic crust obtained using the two approaches differ significantly, not only in velocity magnitude but also structurally. We argue that the results obtained from the acoustic modelling are geologically less plausible and suggest the elastic model as a more reliable representation of the upper-oceanic crust.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉We determine 〈span〉m〈/span〉〈sub〉B〈/sub〉, the original body wave magnitude developed by Gutenberg and Richter over the period 1942–1956, for about 3300 〈span〉M〈/span〉〈sub〉w〈/sub〉 ≥ 6 earthquakes for the period 1988–present using modern broad-band seismograms. The main objective is to extend the database of energy-related parameters by combining 〈span〉m〈/span〉〈sub〉B〈/sub〉 databases for recent and old events. The radiated energy 〈span〉E〈/span〉〈sub〉R_B〈/sub〉 (in erg) computed from 〈span〉m〈/span〉〈sub〉B〈/sub〉 using the Gutenberg & Richter relation $\log {E_{\mathrm{ R}\_\mathrm{ B}}} = 2.4{m_\mathrm{ B}} + 5.8$ agrees very well with 〈span〉E〈/span〉〈sub〉R〈/sub〉 estimated with modern techniques, especially for large deep earthquakes. Thus, 〈span〉E〈/span〉〈sub〉R_B〈/sub〉 is useful as a proxy for 〈span〉E〈/span〉〈sub〉R〈/sub〉 to investigate the global diversity of earthquake characteristics and physics over an extended period of time.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Locating microseismic events is essential for many areas of seismology including volcano and earthquake monitoring and reservoir engineering. Due to the large number of microseismic events in these settings, an automated seismic location method is required to perform real time seismic monitoring. The measurement environment requires a precise and noise-resistant event location method for seismic monitoring. In this paper, we apply Multichannel Coherency Migration (MCM) to automatically locate microseismic events of induced and volcano-tectonic seismicity using sparse and irregular monitoring arrays. Compared to other migration-based methods, in spite of the often sparse and irregular distribution of the monitoring arrays, the MCM can show better location performance and obtain more consistent location results with the catalogue obtained by manual picking. Our MCM method successfully locates many triggered volcano-tectonic events with local magnitude smaller than 0, which demonstrates its applicability on locating very small earthquakes. Our synthetic event location example at a carbon capture and storage site shows that continuous and coherent drilling noise in industrial settings will pose great challenges for source imaging. However, automatic quality control techniques including filtering in the frequency domain and weighting are used to automatically select high-quality data, and can thus effectively reduce the effects of continuous drilling noise and improve source imaging quality. The location performance of the MCM method for synthetic and real microseismic data sets demonstrates that the MCM method can perform as a reliable and automatic seismic waveform analysis tool to locate microseismic events.〈/span〉

Description: 〈span〉〈div〉Abstract〈/div〉In order to investigate the subduction and continental collision of the North China and Yangtze blocks, two magnetotelluric profiles were obtained across the Dabie Orogenic Belt, the Lower Yangtze Depression, and the Jiangnan Orogenic Belt in the central section (Susong-Anqing section) of the middle-lower Reaches of the Yangtze River in China. After data processing and inversions, we obtained electrical models of the crust and upper mantle. The prominent feature revealed by the inversion is an extensive, arched conductive layer that extends from the middle-lower crust to the upper mantle. To the southeast of this layer, another wedge-shaped conductor is located beneath the Lower Yangtze Depression and the Jiangnan Orogenic Belt. In addition, several resistors, which are distributed from the lower crust to the upper mantle, were also revealed by these lines. These resistors are beneath the conductive layer and separated by vertical conductive bands. Based on these electrical structures, we identified several major faults, including the Tanlu Fault in the eastern part of the Dabie Orogenic Belt, the Yangtze Deep Fault, and the main thrust fault in the Lower Yangtze Depression, which are middle-upper crustal faults. In addition, a ‘Crocodile’ structure was revealed by the major faults in the depression, which are connected by a middle-lower crustal detachment and the surrounding resistive strata. Based on the different electrical structures of the three belts and the results of previous studies, we conclude that subsequent to the slab subducting toward the North China Block during the Triassic, the subduction-collision process that occurred in the study area can be divided into three stages. In the first stage, the weak layer and the Yangtze Fold Belt formed during subduction of the Yangtze Block beneath the North China Block, and the Dabie Orogenic Belt formed during the collision process. In the second stage, the slab buckled as deep material upwelled, and the ‘Crocodile’ structure formed due to the weak layer. In the third stage, the breakup and sinking of the slab caused the Yangtze Fold Belt to subside.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉It is well known that for a seismic source on or near a material discontinuity, the moment tensor is ambiguous in the sense that its value depends on the exact location of the source. In this short paper we demonstrate that this ambiguity can be resolved by considering components of the source moment and potency tensors separately. For a general source and any location ‘at’ an interface, certain components of its moment tensor and certain components of its potency tensor can be determined unambiguously. The other components depend on the exact location of the source. To investigate these results it is convenient to allow for an asymmetric moment tensor, where the antisymmetric part is physically equivalent to a torque source. These results are illustrated using ray theory for a source at a strong material contrast.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Earth's large-scale topography and lower mantle structure are linked to past tectonic motions and mantle flow, making it possible to gain insights in the properties of the solid Earth from time-dependent global convection models driven by tectonic reconstructions. Recent work suggests that the amplitude of residual topography, obtained by subtracting models of isostatic topography from total topography, may be up to ∼1,000 m at spherical harmonic degree two (wavelength ∼13,300 km) and 〉 200 m at spherical harmonic degree 12 (wavelength ∼3,000 km). The amplitude of dynamic topography and the structure of the lower mantle predicted by time-dependent forward mantle flow models both depend on the physical assumption and on model parameters. Here, I investigate the consequence of using the Boussinesq, extended Boussinesq or truncated anelastic liquid approximation (TALA) in time-dependent flow models for predicting present-day mantle structure and dynamic topography; I characterise the sensitivity of the spectral distribution of dynamic topography amplitude to the boundary conditions and model setup for the computation of dynamic topography; and I investigate the sensitivity of model results to parameters including the depth- and temperature-dependence of viscosity, the model initial age, and the density of the basal layer. Extended Boussinesq and TALA models are preferred to Boussinesq models that overpredict the volume of lower mantle slabs and the amplitude of long-wavelength dynamic topography. The correlation between dynamic topography and residual topography models at degrees one-three generally ranges between 0.4–0.5. The flow models better predict the geographical location of Large Low Shear Velocity Provinces (LLSVPs) than that of lower mantle slabs, which cover smaller areas in map view. I show that preserving shallow lateral viscosity variations in the computation of dynamic topography increases the amplitude of dynamic topography for wavelengths 〉 6,000 km. Parameter trade-offs exist to fit both deep and surface constraints. The best-fit model cases considered either a moderately dense basal layer (approximately 2 per cent denser than ambient mantle) or weak temperature- and pressure-dependence of lower mantle viscosity. The fit to present-day constraints does not significantly deteriorate when extending the reconstruction of mantle flow from 200 Ma to 410 Ma, reflecting that seismic tomography models capture the history of mantle convection over the last 200–250 Myr, and suggesting that paleogeographically-constrained mantle flow models should be compared to time-dependent surface geological constraints.〈/span〉

Description: 〈span〉〈div〉SUMMARY〈/div〉Guided waves in a water layer overlaying an elastic half-space are known as normal modes. They are often present in seismic recordings at long offsets in shallow-water environment and generally considered coherent noise. The normal modes, however, carry important information about the near-surface and, as demonstrated by a number of authors, can be used to obtain the shallow velocity model. There is a growing evidence that the latter needs not to be isotropic due to various geological reasons. Motivated by that, we consider the normal-mode propagation in case the elastic half-space exhibits orthorhombic anisotropy. We derive the period equation that describes the normal-mode phase velocity dispersion. To simplify the complicated expression, we present acoustic and ellipsoidal orthorhombic approximations. We also outline the approach towards the group velocity and group azimuth calculation and apply it to the ellipsoidal case to obtain concise and intuitive expressions. Using numerical test, we study the relation between phase and group domains in elastic orthorhombic case. The deviation between velocities and azimuths in these domains is the strongest for low frequencies and it rapidly decreases with increasing frequency. For higher frequencies, the anisotropy effects of the underlaying half-space are barely detectable since the observed signal is composed mainly of the direct acoustic wave, resulting in the two domains being nearly indistinguishable.〈/span〉