#### Results in Journal Geophysical Journal International: 17,119

##### (searched for: journal_id:(1130952))
Page of 343
Articles per Page
by
Show export options
Select all
Mark David Wigh, , Arne Døssing
Published: 13 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab371

Abstract:
Summary We investigate if it is theoretically possible to discriminate between unexploded ordnance (UXO) and non-UXO sources by modelling the magnetic dipole moment for ferrous objects of different shapes and sizes. This is carried out by approximating the volumetric demagnetization factors of rectangular prisms, representing shapes similar to a long rod or flat steel plate. By modelling different UXO as prolate spheroids the demagnetization factors can be determined which can be compared with the magnetic response of a prism. The inversion is carried out in a probabilistic framework, where the UXO forward model and the non-UXO forward model are assigned individual prior models in terms of shape, size, orientation and remanent magnetization of the object. 95 independent realizations of the prism prior model are generated to make 95 synthetic anomalies exemplifying non-UXO objects, which are inverted for using the UXO model. It is investigated if an identical magnetic moment can be produced between the two models and how well resolved the magnetic moment is in terms of the measured anomaly. The case study is carried out in two steps where we first have little prior information of expected UXO properties and another where a UXO prior is introduced with expected values of aspect ratio and size of 24 different UXO, that are often encountered in the North Sea. With no prior information of expected UXO, discrimination is at many times implausible, unless elongated rod prism objects are considered, where the magnetic moment often can not be reproduced by a spheroid. Introducing the UXO prior we achieve a much better discrimination rate when using the list of expected UXO properties. By using the UXO prior we can account for a much higher remanent magnetization allowed in the prior, and still achieve high discrimination capabilities in comparison to a case with no UXO prior.
A Michel, J-P Boy
Published: 13 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab369

Abstract:
Summary Long term deformations strongly depend on the Earth model and its rheological parameters, and in particular its viscosity. We give the general theory and the numerical scheme to compute them for any spherically non rotating isotropic Earth model with linear rheology, either elastic or viscoelastic. Although the Laplace transform is classically used to compute viscoelastic deformation, we choose here instead, to implement the integration with the Fourier transform in order to take advantage of the Fast Fourier Transform algorithm and avoid some of the Laplace transform mathematical difficulties. We describe the methodology to calculate deformations induced by several geophysical signals regardless of whether they are periodic or not, especially by choosing an adapted time sampling for the Fourier transform. As examples, we investigate the sensitivity of the displacements due to long period solid Earth tides, Glacial Isostatic Adjustment (GIA), and present-day ice melting, to anelastic parameters of the mantle. We find that the effects of anelasticity are important for long period deformation and relatively low values of viscosities for both Maxwell and Burgers models. We show that slight modifications in the rheological models could significantly change the amplitude of deformation but also affect the spatial and temporal pattern of the signal to a lesser extent. Especially, we highlight the importance of the mantle anelasticity in the low degrees deformation due to present-day ice melting and encourage its inclusion in future models.
Published: 8 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab366

Abstract:
Stagnant-lid convection, where subduction and surface plate motion is absent, is common among the rocky planets and moons in our solar system, and likely among rocky exoplanets as well. How stagnant-lid planets thermally evolve is an important issue, dictating not just their interior evolution but also the evolution of their atmospheres via volcanic degassing. On stagnant-lid planets, the crust is not recycled by subduction and can potentially grow thick enough to significantly impact convection beneath the stagnant lid. We perform numerical models of stagnant-lid convection to determine new scaling laws for convective heat flux that specifically account for the presence of a buoyant crustal layer. We systematically vary the crustal layer thickness, crustal layer density, Rayleigh number, and Frank-Kamenetskii parameter for viscosity to map out system behavior and determine the new scaling laws. We find two end-member regimes of behavior: a “thin crust limit,” where convection is largely unaffected by the presence of the crust, and the thickness of the lithosphere is approximately the same as it would be if the crust were absent; and a “thick crust limit,” where the crustal thickness itself determines the lithospheric thickness and heat flux. Scaling laws for both limits are developed and fit the numerical model results well. Applying these scaling laws to rocky stagnant-lid planets, we find that the crustal thickness needed for convection to enter the thick crust limit decreases with increasing mantle temperature and decreasing mantle reference viscosity. Moreover, if crustal thickness is limited by the formation of dense eclogite, and foundering of this dense lower crust, then smaller planets are more likely to enter the thick crust limit because their crusts can grow thicker before reaching the pressure where eclogite forms. When convection is in the thick crust limit, mantle heat flux is suppressed. As a result, mantle temperatures can be elevated by 100s of degrees K for up to a few Gyrs in comparison to a planet with a thin crust. Whether convection enters the thick crust limit during a planet’s thermal evolution also depends on the initial mantle temperature, so a thick, buoyant crust additionally acts to preserve the influence of initial conditions on stagnant-lid planets for far longer than previous thermal evolution models, which ignore the effects of a thick crust, have found.
, Min Chen, Jieyuan Ning, Tiezhao Bao, Ross Maguire, Megan P Flanagan, Tong Zhou
Published: 4 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab361

Abstract:
Summary The detailed structure near the 410-km discontinuity provides key constraints of the dynamic interactions between the upper mantle and the lower mantle through the mantle transition zone (MTZ) via mass and heat exchange. Meanwhile, the temperature of the subducting slab, which can be derived from its fast wave speed perturbation, is critical for understanding the mantle dynamics in subduction zones where the slab enters the MTZ. Multipathing, i.e. triplicated, body waves that bottom near the MTZ carry rich information of the 410-km discontinuity structure and can be used to constrain the discontinuity depth and radial variations of wave speeds across it. In this study, we systematically analyzed the tradeoff between model parameters in triplication studies using synthetic examples. Specifically, we illustrated the necessity of using array normalized amplitude. Two 1-D depth profiles of the wave speed below the Tatar Strait of Russia in the Kuril subduction zone are obtained. We have observed triplications due to both the 410-km discontinuity and the slab upper surface. And, seismic structures for these two interfaces are simultaneously inverted. Our derived 410-km discontinuity depths for the northern and southern regions are at 420$\pm$15km and 425$\pm$15 km, respectively, with no observable uplift. The slab upper surface is inverted to be located about 50–70 km below the 410-km discontinuity. This location is between the depths of the 1 per cent-2 per cent P-wave speed contours of a regional 3-D FWI model, but we found twice the wave speed perturbation amplitude. A wave speed increase of 3.9 per cent-4.6 per cent within the slab, compared to 2.0 per cent-2.4 per cent from the 3-D full-waveform inversion model, is necessary to fit the waveforms with the shortest period of 2 seconds, indicating that high-frequency waves are required to accurately resolve the detailed structures near the MTZ.
, Jiří Zahradník, Efthimios Sokos, František Gallovič
Published: 3 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab359

Abstract:
Summary A dynamic finite-fault source inversion for stress and frictional parameters of the Mw 6.3 2017 Lesvos earthquake is carried out. The mainshock occurred on June 12, offshore the southeastern coast of the Greek island of Lesvos in the north Aegean Sea. It caused 1 fatality, 15 injuries, and extensive damage to the southern part of the island. Dynamic rupture evolution is modeled on an elliptic patch, using the linear slip-weakening friction law. The inversion is posed as a Bayesian problem and the Parallel Tempering Markov Chain Monte Carlo algorithm is used to obtain posterior probability distributions by updating the prior distribution with progressively more constraints. To calculate the first posterior distribution, only the constraint that the model should expand beyond the nucleation patch is used. Then, we add the constraint that the model should reach a moment magnitude similar to that obtained from our centroid moment tensor inversion. For the final posterior distribution, 15 acceleration records from Greek and Turkish strong motion networks at near regional distances ($\approx 30 - 150$ km) in the frequency range of 0.05–0.15 Hz are used. The three posterior distributions are compared to understand how much each constraint contributes to resolving different quantities. The most probable values and uncertainties of individual parameters are also calculated, along with their mutual trade-offs. The features best determined by seismograms in the final posterior distribution include the position of the nucleation region, the mean direction of rupture (towards WNW), the mean rupture speed (with 68 per cent of the distribution lying between 1.4–2.6 km/s), radiated energy (12–65 TJ), radiation efficiency (0.09–0.38), and the mean stress drop (2.2–6.5 MPa).
I Cho, K Yoshida, H Uebayashi
Published: 3 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab358

Abstract:
Summary The applicability of rotational seismology to the general wavefield of microtremors is theoretically demonstrated based on a random process model of a two-dimensional wavefield. We show the effectiveness of taking the rotations (i.e. spatial differentiation) of microtremor waveforms in separating the Rayleigh and Love waves in a wavefield where waves are simultaneously arriving from various directions with different intensities. This means that a method based on rotational seismology (a rotational method) is capable of separating Rayleigh and Love waves without adopting a specific array geometry or imposing a specific assumption on the microtremor wavefield. This is an important feature of a rotational method because the spatial autocorrelation method, a conventional approach for determining phase velocities in microtremor array surveys, requires either the use of a circular array or the assumption of an isotropic wavefield (i.e. azimuthal averaging of correlations is required). Derivatives of the spatial autocorrelation method additionally require the assumption that Rayleigh and Love waves are uncorrelated. We also show that it is possible to apply a rotational method to determine the characteristics of Love waves based on a simple three-point microtremor array that consists of translational (i.e. ordinary) three-component sensors. In later sections, we assume realistic data processing for microtremor arrays with translational sensors to construct a theoretical model to evaluate the effects of approximating spatial differentiation via finite differencing (i.e. array-derived rotation, ADR) and the effects of incoherent noise on analysis results. Using this model, it is shown that in a short-wavelength range compared to the distance for finite differencing (e.g. $\lambda < 3h$, where $\lambda$ and h are the wavelength and distance for finite differencing, respectively), the leakage of unwanted wave components can determine the analysis limit. It is also shown that in a long-wavelength range (e.g. $\lambda > 3h$), the signal intensity gradually decreases, and thus the effects of incoherent noise increase (i.e. the signal-to-noise ratio decreases) and determine the analysis limit. We derive the relation between the signal-to-noise ratio and wavelength. Although the analysis results quantitatively depend on the array geometry used for finite differencing, the qualitative understanding supported by mathematical expressions with a physically clear meaning can serve as a guideline for the treatment of data obtained from ADR.
W Zürn, , R Widmer-Schnidrig, P Duffner,
Published: 2 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab336

Abstract:
Summary Tilting of the ground due to loading by the variable atmosphere is known to corrupt very long-period horizontal seismic records (below 10 mHz) even at the quietest stations. At BFO (Black Forest Observatory, SW-Germany) the opportunity arose to study these disturbances on a variety of simultaneously operated state-of-the-art broadband sensors. A series of time windows with clear atmospherically caused effects was selected and attempts were made to model these “signals” in a deterministic way. This was done by simultaneously least squares fitting the locally recorded barometric pressure and its Hilbert transform to the ground accelerations in a bandpass between 100 and 3600 s periods. Variance reductions of up to 97 per cent were obtained. We show our results by combining the “specific pressure induced accelerations” for the two horizontal components of the same sensor as vectors on a horizontal plane, one for direct pressure and one for its Hilbert transform. It turned out that at BFO the direct pressure effects are large, strongly position dependent, and largely independent of atmospheric events for instruments installed on piers, while three posthole sensors are only slightly affected. The infamous “cavity effects” are invoked to be responsible for these large effects on the pier sensors. On the other hand, in the majority of cases all sensors showed very similar magnitudes and directions for the vectors obtained for the regression with the Hilbert transform, but highly variable from event to event especially in direction. Therefore this direction most certainly has to do with the gradient of the pressure field moving over the station which causes a larger scale deformation of the crust. The observations are very consistent with these two fundamental mechanisms of how fluctuations of atmospheric surface pressure causes tilt noise. The results provide a sound basis for further improvements of the models for these mechanisms. The methods used here can already help to reduce atmospherically induced noise in long period horizontal seismic records .
Published: 2 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab352

Abstract:
Summary Palaeomagnetic records from geological archives provide significant information about the nature of geomagnetic polarity reversals; however, there are few detailed palaeomagnetic records of pre-Pleistocene reversals. The lower Mammoth Subchron boundary (late Pliocene) is recorded in a 10-m interval of a marine succession deposited at high accumulation rates (9–66 cm/kyr) in the Boso Peninsula, central Japan. Here, we report a continuous palaeomagnetic record of the lower, normal to reverse boundary interval of the Mammoth Subchron, including the geomagnetic field direction and relative palaeointensity, with an average temporal resolution of ca 800 years. A hybrid method of thermal demagnetization at 200° C and progressive alternating field demagnetization were used to effectively extract the primary palaeomagnetic component, which is carried by magnetite. The lower Mammoth transition is characterized by palaeomagnetic direction of instability and decay of the relative palaeointensity, and occurred from late Marine Isotope Stage MG3 (3351 ka) to MG2 (3336 ka) or MG1 (3331 ka), spanning 15–20 kyr. Virtual geomagnetic poles (VGPs), calculated from primary palaeomagnetic directions, rapidly rebounded twice from southern latitudes to northern latitudes within the transition. In contrast to the complex lower Mammoth reversal behavior recorded in the Boso Peninsula succession, records from a lava sequence in O'ahu (Hawai'i) reveal a rebound following a 180° directional change, and those from a marl succession in Sicily (Italy) indicate a single rapid directional change. Diverse geomagnetic field evolution among these three sections is reflected resolution difference among the records likely in combination with an influence of non-axial dipole field.
Published: 1 September 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab355

Abstract:
Summary This paper presents a new sparse inversion method based on L1 norm regularization for 3D magnetic data. In isolation, L1 norm regularization yields model elements which are unconstrained by the input data to be exactly zero, leading to a sparse model with compact and focused structure. Here, we complement the L1 norm with a penalty minimizing total variation, the L1 norm of the model gradients; it is expected that the sharp boundaries of the subsurface structure are not compromised by incorporating this penalty. Although this penalty is widely used in the geophysical inversion studies, it is often replaced by an alternative quadratic penalty to ease solution of the penalized inversion problem; in this study, the original definition of the total variation, i.e., form of the L1 norm of the model gradients, is used. To solve the problem with this combined penalty of L1 norm and total variation, this study introduces alternative direction method of multipliers (ADMM), which is a primal-dual optimization algorithm that solves convex penalized problems based on the optimization of an augmented Lagrange function. To improve the computational efficiency of the algorithm to make this method applicable to large-scale magnetic inverse problems, this study applies matrix compression using the wavelet transform and the preconditioned conjugate gradient method. The inversion method is applied to both synthetic tests and real data, the synthetic tests demonstrate that, when subsurface structure is blocky, it can be reproduced almost perfectly.
Yishan Song, Lian-Feng Zhao, Xiao-Bi Xie, Xiao Ma, Guilin Du, Xiaofeng Tian, Zhen-Xing Yao
Published: 30 August 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab356

Abstract:
Summary On 2019 March 21, an explosion accidentally occurred at a chemical plant in Xiangshui, Yancheng city, Jiangsu province, China. Using broad-band digital seismic data from East China, South Korea and Japan, we investigate properties of the Xiangshui explosion as well as 2 nearby chemical explosions and 4 nearby natural earthquakes in Jiangsu Province, East China. From Lg and Rayleigh waves recorded by regional networks, both body-wave magnitude mb(Lg) and surface-wave magnitude Ms(Rayleigh) are calculated for these events. The magnitudes of the Xiangshui explosion are mb(Lg) = 3.39 ± 0.24 and Ms = 1.95 ± 0.27, respectively. Both the empirical magnitude-yield relation for buried explosion and empirical yield-crater dimension relation for open pit explosion are adopted for investigating the explosive yield. The result from the yield-crater dimension relation is approximately 492 ton, which is consistent with the ground truth and considerably larger than that from the buried source model. This also reveals that, for Xiangshui explosion, the explosion to seismic energy conversion rate is approximately one third compared to a similar sized fully confined explosion. By comparing the body-wave and surface-wave magnitudes from explosions and nearby earthquakes, we find that the mb:Ms discriminant calculated at regional distances cannot properly distinguish explosions from natural earthquakes. However, the P/S spectral ratios Pg/Lg, Pn/Lg and Pn/Sn from the same dataset can be good discriminants for identifying explosions from earthquakes.
, Jean Letort, Fabrice Cotton, Graeme Weatherill, Matthieu Sylvander, Soumaya Latour
Published: 30 August 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab348

Abstract:
Summary Well-constrained earthquake depth estimations are important for seismic hazard determination. As local networks of the East-African Rift are usually too sparse for reliable depth estimations, we used detections of pP and sP phase arrivals (the so-called depth phases) at teleseismic distance to constrain earthquake depths in this region. We rely on a fully automatic Cepstral analysis approach (Letort et al., 2015), first validated at the global scale using the ISC-EHB catalogue, then applied on the East-African seismicity. We investigated 9575 earthquakes from magnitude 2 since 2005 which allows us to constrain the depth estimation of 584 events with magnitude mainly above 3.5, complemented by 139 reliable depth estimations from previous studies based on teleseismic data as well. To ensure a final catalogue as complete as possible, we also identified from regional catalogues 113 earthquakes assumed to be well constrained, based on network geometry empirical criteria. Thanks to this study, we finally propose new earthquake depth distributions for the seismic source zonation defined by Poggi et al. (2017), in order to estimate the seismic hazard of the East African Rift region. Including those new distributions in the source models leads to significant changes of seismic hazard assessments results.
Erratum
, M Ravi Kumar, Santosh Kumar
Published: 30 August 2021
Geophysical Journal International, Volume 227, pp 1937-1937; https://doi.org/10.1093/gji/ggab335

Abstract:
Geophysical Journal International, Volume 226, Issue 3, September 2021, Pages 1800–1813, https://doi.org/10.1093/gji/ggab188.
Nicolás Castro-Perdomo, Renier Viltres, Frédéric Masson, Yann Klinger, Shaozhuo Liu, Maher Dhahry, Patrice Ulrich, Jean-Daniel Bernard, Rémi Matrau, Abdulaziz Alothman, et al.
Published: 28 August 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab353

Abstract:
Summary Although the Dead Sea Transform fault system has been extensively studied in the past, little has been known about the present-day kinematics of its southernmost portion that is offshore in the Gulf of Aqaba. Here we present a new GPS velocity field based on three surveys conducted between 2015 and 2019 at 30 campaign sites, complemented by 11 permanent stations operating near the gulf coast. Interseismic models of strain accumulation indicate a slip rate of $4.9^{+0.9}_{-0.6}~mm/yr$ and a locking depth of $6.8^{+3.5}_{-3.1}~km$ in the gulf’s northern region. Our results further indicate an apparent reduction of the locking depth from the inland portion of the Dead Sea Transform towards its southern junction with the Red Sea rift. Our modelling results reveal a small systematic left-lateral residual motion that we postulate is caused by, at least in part, late postseismic transient motion from the 1995 MW7.2 Nuweiba earthquake. Estimates of the moment accumulation rate on the main faults in the gulf, other than the one that ruptured in 1995, suggest that they might be near the end of their current interseismic period, implying elevated seismic hazard in the gulf area.
, , H Steeb
Published: 28 August 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab354

Abstract:
Summary In this work, we propose a hydro-mechanical simulation model to study the strong interaction of fluid flow and fracture deformation under in-situ stress conditions. The general model is reduced under physics-based assumptions to provide an efficient numerical approach for inverse analysis of experimental studies and is applied to experimental field data obtained from hydraulic tests conducted at the Grimsel Test Site (GTS), Switzerland. The present set of hydro-mechanical measurement data provides not only valuable information about the transient pressure and flow evolution but also the transient change of fracture deformation. We aim to introduce a strongly-coupled hydro-mechanical model to numerically characterize the fractured reservoir based on experimental data below the limit of hydraulically induced irreversible changes of the reservoir’s properties. Insights into the leading mechanisms of flow processes throughout hydraulic testing under in-situ conditions are then gained by best numerical fits of the measurement data. Based on the experimental and numerical findings, this study emphasizes the importance of a consistent consideration of local and non-local fracture deformation throughout inverse analysis of hydraulic testing data to a) better understand hydro-mechanical flow processes in fractured reservoirs and b), to increase the estimation quality of hydraulic properties of tested fractures.
, Jim Mori
Published: 28 August 2021
Geophysical Journal International, Volume 228, pp 387-395; https://doi.org/10.1093/gji/ggab349

Abstract:
SUMMARY Detecting P-wave onsets for online processing is an important component for real-time seismology. As earthquake early warning systems around the world come into operation, the importance of reliable P-wave detection has increased, since the accuracy of the earthquake information depends primarily on the quality of the detection. In addition to the accuracy of arrival time determination, the robustness in the presence of noise and the speed of detection are important factors in the methods used for the earthquake early warning. In this paper, we tried to improve the P-wave detection method designed for real-time processing of continuous waveforms. We used the new Tpd method, and proposed a refinement algorithm to determine the P-wave arrival time. Applying the refinement process substantially decreases the errors of the P-wave arrival time. Using 606 strong motion records of the 2011 Tohoku earthquake sequence to test the refinement methods, the median of the error was decreased from 0.15 to 0.04 s. Only three P-wave arrivals were missed by the best threshold. Our results show that the Tpd method provides better accuracy for estimating the P-wave arrival time compared to the STA/LTA method. The Tpd method also shows better performance in detecting the P-wave arrivals of the target earthquakes in the presence of noise and coda of previous earthquakes. The Tpd method can be computed quickly, so it would be suitable for the implementation in earthquake early warning systems.
, , Laura Behling
Published: 27 August 2021
Geophysical Journal International, Volume 228, pp 275-290; https://doi.org/10.1093/gji/ggab333

Abstract:
SUMMARY Induced polarization (IP) is an acknowledged method in ore exploration and can be applied to evaluate the metal content in dumps containing the residues of ore processing facilities. Existing models explain the relationships between ore content and grain size of the ore particles with IP parameters. However, the models assume full water saturation of the ore containing samples, which is often not the case in field conditions at dump sites. Hence, our study investigates the effect of desaturation on the resulting IP signal. We used six different sand–pyrite mixtures with varying amount and grain sizes of the pyrite particles. Evaporative drying desaturated the samples. Complex conductivity spectra were recorded in the frequency range between 0.02 and 1000 Hz at certain saturation levels. The resulting spectra indicate an decrease of the conductivity amplitude with progressing desaturation. This effect agrees with the second empirical Archie equation. The saturation exponent of the conductivity amplitude shows values slightly larger than one. The measured spectra were processed by a Debye decomposition. We observe a nearly constant total chargeability during desaturation. This finding is in agreement with existing models that relate the total chargeability to the metal content in the sample. However, the mean relaxation time decreases remarkably during the drying process, whereas the normalized relaxation time, which considers the ratio between the mean relaxation time and the resistivity of the embedding material, does not indicate any dependence on water saturation. This behaviour contradicts existing models that predict a decreasing relaxation time with increasing water salinity, which results from evaporative drying. In order to explain the observed effects, we propose a conceptional model that compares a mixture of pyrite particles in an embedding material (sand, water and air) with an electrical RC circuit. The pyrite grains behave as small condensers that are charged and discharged via the conductive background material. According to this simple physical model, the relaxation time is proportional to the resistivity of the embedding material. A resistivity increase while desaturation causes an increase of relaxation time as observed in our experiments. This conceptional model is in good agreement with other experiments that change the resistivity in the background material by varying water salinity or clay content. The capacitive behaviour of the dispersed particles is characterized by the normalized relaxation time that depends on the grain size and the volume content of the pyrite particles.
Jian Wen, Jiankuan Xu, Xiaofei Chen
Published: 27 August 2021
Geophysical Journal International, Volume 228, pp 134-146; https://doi.org/10.1093/gji/ggab346

Abstract:
SUMMARY The stress drop is an important dynamic source parameter for understanding the physics of source processes. The estimation of stress drops for moderate and small earthquakes is based on measurements of the corner frequency ${f_c}$, the seismic moment ${M_0}$ and a specific theoretical model of rupture behaviour. To date, several theoretical rupture models have been used. However, different models cause considerable differences in the estimated stress drop, even in an idealized scenario of circular earthquake rupture. Moreover, most of these models are either kinematic or quasi-dynamic models. Compared with previous models, we use the boundary integral equation method to simulate spontaneous dynamic rupture in a homogeneous elastic full space and then investigate the relations between the corner frequency, seismic moment and source dynamic parameters. Spontaneous ruptures include two states: runaway ruptures, in which the rupture does not stop without a barrier, and self-arresting ruptures, in which the rupture can stop itself after nucleation. The scaling relationships between ${f_c}$, ${M_0}$ and the dynamic parameters for runaway ruptures are different from those for self-arresting ruptures. There are obvious boundaries in those scaling relations that distinguish runaway ruptures from self-arresting ruptures. Because the stress drop varies during the rupture and the rupture shape is not circular, Eshelby's analytical solution may be inaccurate for spontaneous dynamic ruptures. For runaway ruptures, the relations between the corner frequency and dynamic parameters coincide with those in the previous kinematic or quasi-dynamic models. For self-arresting ruptures, the scaling relationships are opposite to those for runaway ruptures. Moreover, the relation between ${f_c}$ and ${M_0}$ for a spontaneous dynamic rupture depends on three factors: the dynamic rupture state, the background stress and the nucleation zone size. The scaling between ${f_c}$ and ${M_0}$ is ${f_c} \propto {M_0^{ - n}}$, where n is larger than 0. Earthquakes with the same dimensionless dynamic parameters but different nucleation zone sizes are self-similar and follow a ${f_c} \propto {M_0^{ - 1/3}}$ scaling law. However, if the nucleation zone size does not change, the relation between ${f_c}$ and ${M_0}$ shows a clear departure from self-similarity due to the rupture state or background stress.
Published: 24 August 2021
Geophysical Journal International, Volume 227, pp 1768-1783; https://doi.org/10.1093/gji/ggab297

Abstract:
SUMMARY The equivalent source technique is a powerful and widely used method for processing gravity and magnetic data. Nevertheless, its major drawback is the large computational cost in terms of processing time and computer memory. We present two techniques for reducing the computational cost of equivalent source processing: block-averaging source locations and the gradient-boosted equivalent source algorithm. Through block-averaging, we reduce the number of source coefficients that must be estimated while retaining the minimum desired resolution in the final processed data. With the gradient-boosting method, we estimate the sources coefficients in small batches along overlapping windows, allowing us to reduce the computer memory requirements arbitrarily to conform to the constraints of the available hardware. We show that the combination of block-averaging and gradient-boosted equivalent sources is capable of producing accurate interpolations through tests against synthetic data. Moreover, we demonstrate the feasibility of our method by gridding a gravity data set covering Australia with over 1.7 million observations using a modest personal computer.
Y X Zhao, Y Li, N Wu
Published: 24 August 2021
Geophysical Journal International, Volume 228, pp 119-133; https://doi.org/10.1093/gji/ggab345

Abstract:
SUMMARY As a data-driven approach, the performance of deep learning models depends largely on the quantity and quality of the training data sets, which greatly limits the application of deep learning to tasks with small data sets. Unfortunately, sometimes we need to use limited small data sets to complete our tasks, such as distributed acoustic sensing (DAS) data denoising. However, using a small data set to train the network may cause overfitting, resulting in poor network generalization. To solve this problem, we propose an approach based on the combination of a generative adversarial network and a deep convolutional neural network. First, we used a small noise data set to train a generative adversarial network to generate synthetic noise samples, and then used these synthetic noise samples to augment the noise data set. Next, we used the augmented noise data set and the signal data set obtained through forward modelling to construct a synthetic training set. Finally, a denoising network based on a convolutional neural network was trained on the constructed synthetic training set. Experimental results show that the augmented data set can effectively improve the denoising performance and generalization ability of the network, and the denoising network trained on the augmented data set can more effectively reduce various kinds of noise in the DAS data.
, J Hinderer, F Masson, D Viville, S Pasquet, , J D Bernard, N Lesparre, M C Pierret
Published: 24 August 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab328

Abstract:
Summary Assessing the spatial and temporal heterogeneity in subsurface water storage has strong societal and environmental implications, as it is key to assess the water availability for the ecosystem and society. This challenge is especially significant in mountainous areas, where the local population totally depends on springwater as a freshwater resource, while water storage dynamics is complex to evaluate because it exhibits spatiotemporal heterogeneities on all scales as a result of the topography. In this study, we compare the water balance of a headwater granitic catchment (CWB) with water storage changes assessed from in-situ continuous gravity monitoring using an iGrav superconducting gravimeter (SGWSC) located at the summit of the catchment. We show that SGWSC and CWB exhibit a similar annual cycle, although they deviate in the months following winter peak flow events. We investigate the reasons for these discrepancies using a tank model adjusted to the SG signal. This shows that during these events, the effective discharge in the SG footprint area is much lower than the catchment streamflow. We attribute this difference in the drainage term to a lower contribution of the upper part of the catchment to the generation of peak flow, compared to the lower part.
, Christine A Powell
Published: 23 August 2021
Geophysical Journal International, Volume 228, pp 291-307; https://doi.org/10.1093/gji/ggab337

Abstract:
SUMMARY Phase velocity and azimuthal anisotropy maps for fundamental mode Rayleigh waves are determined for a portion of the central United States including the seismically active Reelfoot Rift (RFR) and the enigmatic Illinois Basin. Dense seismic array installations of the Northern Embayment Lithosphere Experiment, the EarthScope transportable array and the Ozarks Illinois Indiana Kentucky array allow a detailed investigation of phase velocity and anisotropy in a broad period range (20–100s).We obtain more than 12 000 well-constrained, unique two-station paths from teleseismic events. The two-station method is used to determine dispersion curves and these are inverted for isotropic phase velocity maps and azimuthal anisotropy maps for each period. The presence of fast phase velocities at lower crustal and uppermost mantle depths is found below the RFR, and Ste. Genevieve and Wabash Valley fault zones. At periods of 30s and higher, the RFR is underlain by slow phase velocities and is flanked to the NW and SE by regions of fast velocity. Fast phase velocities are present below the centre of the Illinois Basin in the period range 75–100s. Anisotropy fast axis orientations display complex patterns for each period and do not trend parallel to the direction of absolute plate motion. Anisotropy fast directions are consistently parallel to the trend of the RFR from 50s to higher periods, suggesting the presence of either frozen-in anisotropic fabric or fabric related to material transport from a recently discovered, pronounced low velocity zone below the Mississippi Embayment.
Published: 23 August 2021
Geophysical Journal International, Volume 228, pp 355-365; https://doi.org/10.1093/gji/ggab344

Abstract:
SUMMARY Correlation detectors are now used routinely in seismology to detect occurrences of signals bearing close resemblance to a reference waveform. They facilitate the detection of low-amplitude signals in significant background noise that may elude detection using energy detectors, and they associate a detected signal with a source location. Many seismologists use the fully normalized correlation coefficient C between the template and incoming data to determine a detection. This is in contrast to other fields with a longer tradition for matched filter detection where the theoretically optimal statistic C2 is typical. We perform a systematic comparison between the detection statistics C and C|C|, the latter having the same dynamic range as C2 but differentiating between correlation and anticorrelation. Using a database of short waveform segments, each containing the signal on a 3-component seismometer from one of 51 closely spaced explosions, we attempt to detect P- and S-phase arrivals for all events using short waveform templates from each explosion as reference signals. We present empirical statistics of both C and C|C| traces and demonstrate that C|C| detects confidently a higher proportion of the signals than C without evidently increasing the likelihood of triggering erroneously. We recall from elementary statistics that C2, also called the coefficient of determination, represents the fraction of the variance of one variable which can be explained by another variable. This means that the fraction of a segment of our incoming data that could be explained by our signal template decreases almost linearly with C|C| but diminishes more rapidly as C decreases. In most situations, replacing C with C|C| in operational correlation detectors may improve the detection sensitivity without hurting the performance-gain obtained through network stacking. It may also allow a better comparison between single-template correlation detectors and higher order multiple-template subspace detectors which, by definition, already apply an optimal detection statistic.
, Jingsheng Ma, Xiaoyang Wu, Dorrik Stow
Published: 20 August 2021
Geophysical Journal International, Volume 228, pp 39-65; https://doi.org/10.1093/gji/ggab334

Abstract:
SUMMARY An evaluation of prospective shale gas reservoir intervals in the Bowland Shale is presented using a wireline log data set from the UK's first shale gas exploration well. Accurate identification of such intervals is crucial in determining ideal landing zones for drilling horizontal production wells, but the task is challenging due to the heterogeneous nature of mudrocks. This heterogeneity leads to stratigraphic variations in reservoir quality and mechanical properties, and leads to complex geophysical behaviour, including seismic anisotropy. We generate petrophysical logs such as mineralogy, porosity, and organic content and calibrate these to the results of core studies. If ‘reservoir quality’ is defined by combined cut-offs relating to these parameters, we find that over 100 m of reservoir quality shale is present in the well, located primarily within the upper section. To examine the link between geophysical signature and rock properties, an isotropic rock physics model is developed, using effective medium theories, to recreate the elastic properties of the shale and produce forward-looking templates for subsequent seismic inversion studies. We find that the mineralogical heterogeneity in the shale has a profound impact on modelled elastic properties, obscuring more discrete changes due to porosity, organic content and water saturation and that the best reservoir quality intervals of the shale bear a distinctive response on rock physics cross-plots. Finally, we consider the density of natural fractures in the shale by developing an anisotropic rock physics model to reflect high-angle fractures observed on micro-imagery logs. We invert crack density using shear wave splitting well log data and find a crack density of up to 4 per cent which correlates well with micro-imagery observations. Our work further supports previous authors’ core-based studies in concluding that the Bowland Shale holds good reservoir characteristics, and we propose that there are multiple intervals within the shale that could be targeted with stacked horizontal wells, should those intervals’ mechanical properties also be suitable and there be adequate stress barriers between to restrict vertical hydraulic fracture growth. Finally, our rock physics templates may provide useful tools in interpreting pre-stack seismic data sets in prospective areas of the Bowland Shale and picking the best locations for drilling wells.
, Richard K Bono, Domenico G Meduri, , Samuel Greenwood, Andrew J Biggin
Published: 20 August 2021
Geophysical Journal International, Volume 228, pp 316-336; https://doi.org/10.1093/gji/ggab342

Abstract:
SUMMARY Elucidating the processes in the liquid core that have produced observed palaeointensity changes over the last 3.5 Gyr is crucial for understanding the dynamics and long-term evolution of Earth’s deep interior. We combine numerical geodynamo simulations with theoretical scaling laws to investigate the variation of Earth’s magnetic field strength over geological time. Our approach follows the study of Aubert et al., adapted to include recent advances in numerical simulations, mineral physics and palaeomagnetism. We first compare the field strength within the dynamo region and on the core–mantle boundary (CMB) between a suite of 314 dynamo simulations and two power-based theoretical scaling laws. The scaling laws are both based on a Quasi-Geostropic (QG) force balance at leading order and a Magnetic, Archimedian, and Coriolis (MAC) balance at first order and differ in treating the characteristic length scale of the convection as fixed (QG-MAC-fixed) or determined as part of the solution (QG-MAC-free). When the data set is filtered to retain only simulations with magnetic to kinetic energy ratios greater than at least two we find that the internal field together with the root-mean-square and dipole CMB fields exhibit power-law behaviour that is compatible with both scalings within uncertainties arising from different heating modes and boundary conditions. However, while the extrapolated intensity based on the QG-MAC-free scaling matches Earth’s modern CMB field, the QG-MAC-fixed prediction shoots too high and also significantly overestimates palaeointensities over the last 3.5 Gyr. We combine the QG-MAC-free scaling with outputs from 275 realizations of core–mantle thermal evolution to construct synthetic true dipole moment (TDM) curves spanning the last 3.5 Gyr. Best-fitting TDMs reproduce binned PINT data during the Bruhnes and before inner core nucleation (ICN) within observational uncertainties, but PINT does not contain the predicted strong increase and subsequent high TDMs during the early stages of inner core growth. The best-fitting models are obtained for a present-day CMB heat flow of 11–16 TW, increasing to 17–22 TW at 4 Ga, and predict a minimum TDM at ICN.
, , Alexis Le Pichon, Junghyun Park, Stephen Arrowsmith, Brian Stump
Published: 20 August 2021
Geophysical Journal International, Volume 228, pp 308-315; https://doi.org/10.1093/gji/ggab338

Abstract:
SUMMARY North Korea conducted its sixth underground nuclear explosion test (${m_\mathrm{ b}}$ 6.3) on 2017 September 3. The underground explosion produced substantial low-frequency atmospheric waves, which were detected by infrasound arrays located up to a distance of 566 km. These infrasound waves are formed by the conversion of seismic energy to acoustic energy across the lithosphere–atmosphere interface. While infrasound records at regional distances produce estimates of ground motion amplitude over spatially extended regions covering about 26 500 km2, 3-D full seismo-acoustic simulations within the lithosphere and atmosphere provide quantitative information about seismo-acoustic energy partitioning. Our results demonstrate the capability of remote infrasound observations combined with 3-D propagation modelling to further develop discrimination methods for underground sources. These results contribute to enhance the confidence of source identification and characterization in nuclear test monitoring research, which is essential for the enforcement of the Comprehensive Nuclear-Test-Ban Treaty.
, Alexey Stovas, , Xingxing Huang, Lingyi He
Published: 20 August 2021
Geophysical Journal International, Volume 228, pp 102-118; https://doi.org/10.1093/gji/ggab341

Abstract:
SUMMARY The monoclinic medium with a horizontal symmetry plane is gradually being studied for seismic anisotropy characterization. The principle goal of this paper is to investigate the effect of fracture parameters (azimuth angle, density, aspect ratio, scale) on the exact and approximate monoclinic anisotropy parameters. We derive the monoclinic porous media based on the Chapman model which accounts for the wave-induced fluid flow and give the expressions of the Thomsen-style anisotropy parameters (nine orthorhombic anisotropy parameters: VP0, VS0, ϵ1, ϵ2, γ1, γ2, δ1, δ2, δ3, three exact monoclinic parameters: ζ1, ζ2, ζ3 and three approximate monoclinic parameters: $\widetilde{\zeta _{1}}, \widetilde{\zeta _{2}}, \widetilde{\zeta _{3}}$). The dependence of Thomsen-style anisotropy parameters associated with azimuth angle between two fracture sets is analysed. The orthorhombic anisotropy parameters and monoclinic anisotropy parameters have the same period (π) on the azimuth angle between two fracture sets. The exact and approximate monoclinic anisotropy parameters responsible for the rotation of the P-wave NMO ellipse have a similar trend versus the azimuth angle, while those responsible for the rotation of the S1- and S2-wave NMO ellipses have significant discriminations. The influence of fracture density, aspect ratio, and scale on the monoclinic parameters are also analysed. The monoclinic anisotropy parameters responsible for the rotation of the P-wave NMO ellipse decrease with fracture density and aspect ratio increasing from 0 to 0.1, while those responsible for the rotation of S1- and S2-wave NMO ellipses increase with the fracture parameters. The fracture density has a bigger influence on the monoclinic anisotropy parameters than the fracture aspect ratio. When saturated with different fluids (water and CO2), the monoclinic parameters have a similar behaviour versus the azimuth angle between two fracture sets.
, Dirk Klaeschen, Christian Berndt, Simona Pierdominici, Peter Hedin
Published: 20 August 2021
Geophysical Journal International, Volume 228, pp 66-77; https://doi.org/10.1093/gji/ggab339

Abstract:
SUMMARY Strong anisotropy of seismic velocity in the Earth's crust poses serious challenges for seismic imaging. Where in situ seismic properties are not available, the anisotropy can be determined from velocity analysis of surface and borehole seismic profiles. This is well established for dense, long-offset reflection seismic data. However, it is unknown how applicable this approach is for sparse seismic reflection data with low fold and short offsets in anisotropic metamorphic rocks. Here, we show that anisotropy parameters can be determined from a sparse 3-D data set at the COSC-1 borehole site in the Swedish Caledonides and that the results agree well with the seismic anisotropy parameters determined from seismic laboratory measurements on core samples. Applying these anisotropy parameters during 3-D seismic imaging improves the seismic image of the high-amplitude reflections especially in the vicinity of the lower part of the borehole. Strong reflections in the resulting seismic data show good correlation with the borehole-derived lithology. Our results aid the interpretation and extrapolation of the seismic stratigraphy of the Lower Seve Nappe in Jämtland and other parts in the Caledonides.
, , QuanCheng Huang, Nicholas Schmerr
Published: 20 August 2021
Geophysical Journal International, Volume 228, pp 78-101; https://doi.org/10.1093/gji/ggab340

Abstract:
SUMMARY We investigated the likelihood of radial anisotropy in the shallow and deep upper mantle, including the mantle transition zone (MTZ) under the Indian Ocean. Seismic anisotropy can be an indicator of mantle deformation through lattice preferred orientation of anisotropic crystals in the mantle. It has thus the potential to illuminate Earth’s dynamic interior, but previous seismic tomography studies have not achieved consensus on the existence of radial anisotropy below ∼250 km depth. We developed a fully nonlinear transdimensional hierarchical Bayesian Markov Chain Monte Carlo approach to invert fundamental and higher mode surface wave dispersion data and applied it to a subset of a global Love and Rayleigh wave data set. We obtained posterior model parameter distributions for shear wave velocity (VS) and radial anisotropy ξ under the Indian Ocean. These posterior model distributions were used to calculate the probability of having radial anisotropy at different depths. We demonstrated that separate inversions of Love and Rayleigh waves yield models compatible with the results of joint inversions within uncertainties. The obtained pattern of VS anomalies agrees with most previous studies. They display negative anomalies along ridges in the uppermost mantle, but those are stronger than for regularized inversions. The Central Indian Ridge and the Southeastern Indian Ridge present velocity anomalies that extend to ∼200 km depth, whereas the Southwestern Indian Ridge seems to have a shallower origin. Weaker, laterally variable velocity perturbations were found at larger depths. The anisotropy models differ more strongly from regularized inversion results and their uncertainties were rather large. We found that anisotropy models from regularized inversions also depend on the chosen parametrization, which is consistent with the existence of a large model null-space. Apart from a fast horizontally polarized shear wave signal in the top 100 km, likely reflecting the horizontal plate motion due to asthenospheric deformation, no clear relation to surface geology was detected in the anisotropy models. We found that, although the anisotropy model uncertainties are rather large, and lateral variations are present, the data generally prefer at least 1 per cent anisotropy in the MTZ with fast vertically polarized shear waves, within errors. Incorporating group velocity data did not help better constrain deep structure by reducing parameter trade-offs. We also tested the effect of prior constraints on the 410- and 660-km topography and found that the undulations of these discontinuities had little effect on the resulting models in our study region.
Published: 17 August 2021
Geophysical Journal International, Volume 228, pp 17-38; https://doi.org/10.1093/gji/ggab331

Abstract:
SUMMARY A numerical simulation of earthquake cycles at the subduction zone of a plate interface was conducted, using a rate- and state-dependent friction law, to examine aseismic sliding processes propagating into a locked plate interface. In the model, a reverse fault is assumed in a 2-D uniform elastic half-space, and relative plate motion is imposed with a constant velocity. Simulated earthquakes occur repeatedly at a shallower seismogenic plate interface with a velocity-weakening frictional property, while stable sliding occurs in the deeper part, in which a velocity-strengthening frictional property is assumed. During interseismic periods, deep stable sliding causes shear stress to concentrate at the deeper edge of the locked plate interface, and a partial drop in stress occurs, resulting in plate detachment. The resulting detachment front propagates upwards along the seismogenic plate interface until an earthquake occurs, causing aseismic sliding with a sliding velocity significantly lower than the imposed relative plate velocity. The propagation velocity of the detachment front is almost constant in each case, and it is proportional to the relative plate velocity and inversely proportional to the effective normal stress. Episodic events of increased slip velocity occur in the latter half of an interseismic period when the characteristic slip distance is small. Updip propagation of episodic slip is arrested by a low stress barrier, and episodic slips backpropagate downdip. During the final few years of the simulation cycle, the average sliding velocity is approximately inversely proportional to the time to occurrence for a large interplate earthquake.
, Abdolhamid Ansari
Published: 17 August 2021
Geophysical Journal International, Volume 228, pp 1-16; https://doi.org/10.1093/gji/ggab330

Abstract:
SUMMARY Magnetization vector inversion (MVI) has attracted considerable attention in recent years since by this inversion both distribution of the magnitude and direction of the magnetization are obtained; therefore, it is easy to distinguish between the magnetic causative bodies especially when magnetic data are affected by different remanent magnetization. In this research, the compact magnetization vector inversion is presented: a 3-D magnetic modelling is proposed from surface data measurements to obtain compact magnetization distribution. The equations are solved in data-space least squares and the algorithm includes a combination of two weights as depth weighting and compactness weighting in the Cartesian system. The re-weighted compactness weighting matrix handles sparsity constraints imposed on the magnitude of magnetization for varying Lp-norms ($0 \le p \le 2$). The low value of the norm leads to more focused or compact inversion, and using a large value of p obtains a smooth model. The method is validated with two synthetic examples, the first is a cube that has significant remanent magnetization and the second consists of two causative cube bodies with significant different magnetization directions at different depths. The case study is the magnetic data of Galinge iron ore deposit (China) that the apparent susceptibility and magnetization directions are reconstructed. The compact model reveals that the results agree with drilling and geological information.
, Zhinong Wang, Dunshi Wu, Ruiqian Cai, Han Wu
Published: 16 August 2021
Geophysical Journal International, Volume 227, pp 1480-1495; https://doi.org/10.1093/gji/ggab284

Abstract:
SUMMARY Compared with surface waves, guided waves are rarely applied in near-surface investigation. The main reason may lie in the complexity of their dispersion curves. Besides, the study and understanding of guided wave dispersion characteristics are now also inadequate and not deep enough. In this paper, we derived the complete theoretical dispersion curves of P–SV-wave and pure P-wave systems in layered media based on the transmission matrix method and obtained the relative displacement amplitude coefficients at the free surface as a function of frequency and phase velocity for both surface and guided waves. By assigning the value of relative displacement amplitude coefficient to the corresponding point (f,v) on dispersion curve, we got a multi-information diagram called relative amplitude dispersion map (RADM). As a unified description of surface and guided waves, RADM not only shows the velocity–frequency relationship but also represents the polarized energy ratio at the free surface by display colours. The accuracy of RADM was proved by synthetic seismic records, in which RADMs fit well with the corresponding dispersion energy of surface and guided waves. In addition, we designed six models with different Poisson's ratio (PR) and different layer numbers for comparison. It shows that the dispersive vertical-to-horizontal amplitude ratio of guided waves is complex and discontinuous in RADM, which brings great difficulty for mode identification and even affects the subsequent inversion. Tests also show that for high PR layers, the trends of guided P–SV-wave dispersion curves are basically consistent with those of pure P wave. With the decrease of PR, dispersion curves of guided P–SV wave gradually deviate from those of pure P wave. However, RADMs can be greatly consistent with the dispersion energy in either case. This is of great significance for the inversion of near-surface P and S velocities by using dispersion relationships of multimode surface and guided waves.
, Li Zhao, Xiaofei Chen
Published: 16 August 2021
Geophysical Journal International, Volume 228, pp 250-258; https://doi.org/10.1093/gji/ggab329

Abstract:
SUMMARY A novel approach to computing the Earth's normal modes is presented. With the Langer approximation incorporated into the framework of the theory of the generalized reflection/transmission coefficients, our method, when applied to calculate the toroidal modes, produces results that agree extremely well with the exact results from Mineos and successfully overcomes the limitations of the WKBJ analysis. Given any 1-D earth model, regardless of the number of discontinuities, our method can perform the calculations reliably with satisfactory accuracy at high frequencies. The success achieved with the toroidal modes encourages us to tackle in a future study the asymptotic computation of the spheroidal modes, especially those high-frequency trapped modes for which the accuracy of Mineos is demonstrably inadequate.
K Gwirtz, M Morzfeld, W Kuang, A Tangborn
Published: 14 August 2021
Geophysical Journal International, Volume 227, pp 2180-2203; https://doi.org/10.1093/gji/ggab327

Abstract:
SUMMARY Geomagnetic data assimilation merges past and present-day observations of the Earth’s magnetic field with numerical geodynamo models and the results are used to initialize forecasts. We present a new ‘proxy model’ that can be used to test, or rapidly prototype, numerical techniques for geomagnetic data assimilation. The basic idea for constructing a proxy is to capture the conceptual difficulties one encounters when assimilating observations into high-resolution, 3-D geodynamo simulations, but at a much lower computational cost. The framework of using proxy models as ‘gate-keepers’ for numerical methods that could/should be considered for more extensive testing on operational models has proven useful in numerical weather prediction, where advances in data assimilation and, hence, improved forecast skill, are at least in part enabled by the common use of a wide range of proxy models. We also present a large set of systematic data assimilation experiments with the proxy to reveal the importance of localization and inflation in geomagnetic data assimilation.
, Hiroyuki Kumagai,
Published: 14 August 2021
Geophysical Journal International, Volume 227, pp 2137-2155; https://doi.org/10.1093/gji/ggab325

Abstract:
SUMMARY Long-period (LP) seismic events have occurred repeatedly at Galeras volcano, Colombia, during the transition from effusive dome formation to explosive Vulcanian eruptions. Since 1989, two types of LP events have been observed there: one characterized by long-lasting, decaying harmonic oscillations (NLP events) and the other by non-harmonic oscillatory features (BLP events). NLP events are attributed to resonances of a dusty gas-filled crack in the magma plugging the eruptive conduit. Sixteen episodes of NLP events occurred at Galeras during 1992–2010, each characterized by systematic temporal variations in the frequencies and quality factors of NLP events. Our and previous estimates of crack model parameters during three of those NLP episodes indicate that the similar temporal variations in crack geometry and fluid properties can be explained by an increase in the ash content within the crack and a decrease in crack volume. We found that NLP events, associated with low SO2 fluxes, are anticorrelated with BLP events, which are accompanied by high SO2 emissions. From our observations and analytical results, we inferred that BLP events are generated by resonances of open cracks in the uppermost magma plug, corresponding to tuffisite veins, that efficiently transfer volcanic gases. After sufficient degassing and densification, the magma plug effectively seals the conduit. The growing overpressure in the deeper magma is then released through a shear fracture along the conduit margin. The intrusion of deeper, vesiculated magma into the shear fracture depressurizes and fragments the magma, producing a dusty gas and triggering the crack resonances that generate NLP events. Our results thus indicate that the evolution of the properties of the magma plug controls the occurrences of BLP and NLP events at Galeras. Although NLP events do not always precede explosive eruptions, they indicate that an important overpressure is building in the shallow conduit.
, Waldemar Jóźwiak, Krzysztof Nowożyński
Published: 13 August 2021
Geophysical Journal International, Volume 227, pp 1917-1936; https://doi.org/10.1093/gji/ggab322

Abstract:
SUMMARY Long-period magnetotelluric soundings were carried out in the eastern part of the Fore-Sudetic Monocline in central Europe to determine the deep geoelectrical structure. It is an important area in the contact zone of the Palaeozoic terranes of Western Europe and the Proterozoic East European Craton. The research area includes the Dolsk fault zone and the Odra fault zone, which are essential geological boundaries of a regional nature separating crustal blocks of various origins. There were conducted 51 soundings on the quasi-regular mesh 150 km in SE and 225 km in NE. The research region is highly urbanized, and hence some of the data were strongly disturbed. Careful processing of data was required, and sometimes measurements had to be repeated at other locations. The collected data allow constructing 3-D models of the electrical conductivity distribution. A parallel version of the ModEM 3-D inversion code was used for modelling. The information contained in each used transfer function was also examined by performing a separate inversion of these functions. The obtained resistivity distribution models confirm the Dolsk and Odra faults' location as postulated previously by geologists. They also show that these very deep faults are reaching the lower crust. However, they do not clearly state how far the contemporary lowering of the Baltica reaches under cover of younger Palaeozoic sediments.
Published: 13 August 2021
Geophysical Journal International, Volume 227, pp 2156-2179; https://doi.org/10.1093/gji/ggab326

Abstract:
SUMMARY Attenuation mechanism of seismic waves derived with an appropriate scattering model is a mandate for high-frequency earthquake ground-motion modelling. The assumption of uniform half-space is not always realistic and can have significant implications on the results obtained. We tested this assumption for the case of Garhwal Himalaya using a two-layered half-space model. To this end, we analysed the 1999 Chamoli earthquake aftershock sequence in Garhwal Himalaya based on the MLTWA method assuming both the uniform and layered models. The data set consists of 384 seismograms from 72 aftershocks recorded by a temporary seismic network of twelve stations. The uniform model cannot fully describe the energy loss by intrinsic and scattering attenuation, expecting bias in the corresponding attenuation estimates which are constrained by comparing results of numerical simulations with analytical solutions. By comparing the coda envelopes of both the analytical and numerical experiments for a suite of models that include the varying degree of scattering and intrinsic attenuation of crust and mantle (half-space), we could obtain improved attenuation estimates. Monte Carlo method was used to numerically solve the radiative transfer equation to deal with multiple isotropic scattering of seismic waves in 3-D heterogeneous acoustic media. The coda envelopes for the uniform model could reasonably fit those for the layered ones. The envelope fits suggest that the intrinsic and scattering coefficients for the uniform model are overestimated, by a factor that is frequency-dependent, as compared to those for the layered earth model. The factor lies in the range 4.0–7.5 (5.75 ± 1.75) below ∼3 Hz and 2.8–4.0 (3.4 ± 0.6) above 3 Hz for intrinsic, and the range 1.1–2.5 (1.8 ± 0.7) for scattering coefficient at 1–12 Hz. Under the assumption of a layered model, we could further reinterpret the apparent frequency-dependent energy loss as a sum of frequency-dependent mantle leakage and frequency-independent intrinsic attenuation. The results in the assumption of a layered model suggest strong scattering in the upper crust and weak mantle leakage in the frequency range 1.5–24 Hz. The small-scale random heterogeneities responsible for the observed scattering likely represent the fluid-filled fractured-crust, as evidenced from seismic tomography in the source area of the 1999 event.
Zhixiang Chen, , Longtan Shao, Shunqun Li, Lingxia Gao
Published: 13 August 2021
Geophysical Journal International, Volume 228, pp 240-249; https://doi.org/10.1093/gji/ggab319

Abstract:
SUMMARY Latent heat of soil water phase transition is an important parameter affecting the prediction of the temperature field in a frozen soil. With the nonlinear change of unfrozen water content in freezing soil, the latent heat appears nonlinearly in a frozen soil. To evaluate the influence of the estimation of phase quantity and phase change interval on the formation of the temperature field in a frozen soil, the thermal parameters at different temperatures are measured, and a soil freezing model test is carried out. Meanwhile, based on the evolutionary trend of the unfrozen water content, a precise method of transforming the latent heat of phase change into specific heat and a segment processing method of latent heat in a frozen soil are proposed. Also, these methods are used in the ABAQUS to simulate the temperature of the soil freezing model. The sensitivity of the estimation of phase quantity and phase change interval on the formation of the temperature field in a frozen soil is analysed based on the numerical calculation values and model test values. The results show that the estimation accuracy of phase change interval greatly affects the prediction accuracy of the temperature field in the warm temperature stage, and the method of transforming latent heat into apparent specific heat can improve the prediction accuracy of the temperature field.
Yongming Lu,
Published: 13 August 2021
Geophysical Journal International, Volume 227, pp 2121-2136; https://doi.org/10.1093/gji/ggab324

Abstract:
SUMMARY Traveltime computations are an important aspect of seismic data processing applications such as traveltime tomography, migration and seismic source localization. Seismic anisotropy is a widespread feature of the Earth. Solutions to the eikonal equation that account for such anisotropy are needed for high-resolution seismic imaging and inversion. The fast sweeping method (FSM) has been widely used in computing the first-arrival traveltimes for anisotropic media because it does not need to expand the wave front from the point of the smallest traveltime. To apply FSM on strong anisotropic media, one has to solve the slowness equation derived from the Christoffel equation. All the previous developed FSM methods transform the quartic coupled slowness surface equation of quasi-P (qP) and quasi-SV (qSV) waves to the quartic equation in terms of the unknown traveltime, then numerically solve this quartic equation to compute the first-arrival traveltimes of the qP waves. However, the computational cost is significantly increased due the numerically solving the quartic equation, especially for the 3-D problems. In this study, we find a way to transform the quartic slowness equation into a quadratic one if a specific triangular-pyramid stencil around a target point is used. As the quadratic equation has the analytical solution and does not need a numerical solver, the computational efficiency of the scheme is greatly improved. We apply this methodology to develop an efficient 3-D FSM to compute the first-arrival traveltimes for qP waves in 3-D vertical transversely isotropic (VTI) media. We use both layered VTI model and complex VTI model to demonstrate the efficiency of the proposed method to obtain accurate traveltimes in 3-D VTI media involving strong anisotropy.
Published: 13 August 2021
Geophysical Journal International, Volume 227, pp 2095-2120; https://doi.org/10.1093/gji/ggab323

Abstract:
SUMMARY We present an inversion algorithm tailored for point gravity data. As the data are from multiple surveys, it is inconsistent with regards to spacing and accuracy. An algorithm design objective is the exact placement of gravity observations to ensure no interpolation of the data is needed prior to any inversion. This is accommodated by discretization using an unstructured tetrahedral finite-element mesh for both gravity and density with mesh nodes located at all observation points and a first-order system least-squares (FOSLS) formulation for the gravity modelling equations. Regularization follows the Bayesian framework where we use a differential operator approximation of an exponential covariance kernel, avoiding the usual requirement of inverting large dense covariance matrices. Rather than using higher order basis functions with continuous derivatives across element faces, regularization is also implemented with a FOSLS formulation using vector-valued property function (density and its gradient). Minimization of the cost function, comprised of data misfit and regularization, is achieved via a Lagrange multiplier method with the minimum of the gravity FOSLS functional as a constraint. The Lagrange variations are combined into a single equation for the property function and solved using an integral form of the pre-conditioned conjugate gradient method (I-PCG). The diagonal entries of the regularization operator are used as the pre-conditioner to minimize computational costs and memory requirements. Discretization of the differential operators with the finite-element method (FEM) results in matrix systems that are solved with smoothed aggregation algebraic multigrid pre-conditioned conjugate gradient (AMG-PCG). After their initial setup, the AMG-PCG operators and coarse grid solvers are reused in each iteration step, further reducing computation time. The algorithm is tested on data from 23 surveys with a total of 6519 observation points in the Mt Isa–Cloncurry region in north–west Queensland, Australia. The mesh had about 2.5 million vertices and 16.5 million cells. A synthetic case was also tested using the same mesh and error measures for localized concentrations of high and low densities. The inversion results for different parameters are compared to each other as well as to lower order smoothing. Final inversion results are shown with and without depth weighting and compared to previous geological studies for the Mt Isa–Cloncurry region.
Published: 12 August 2021
Geophysical Journal International, Volume 227, pp 1905-1916; https://doi.org/10.1093/gji/ggab321

Abstract:
SUMMARY Simplified solutions of the Bloch equation can lead to inaccurate estimates of hydrogeological parameters from surface nuclear magnetic resonance measurements. Even for single pulse measurements, using simplified forward models is common practice because of the computational intensity of obtaining the full-Bloch solution. These challenges are exacerbated for multipulse sequences. We show parallelizing the full-Bloch solver on a Graphics Processing Unit reduces the solve time by three orders of magnitude. Further optimizations by numerical, analytical and hybrid solutions yield an additional 3× speed up. We simulate the full-Bloch physics for free-induction decay, spin-echo and pseudo-saturation recovery excitation schemes for an unprecedented range of physical scenarios. We explore the time-dependence and relaxation time sensitivity in these solution spaces. Characterizing the solution spaces with polynomials of the relaxation times, the solutions can be rapidly reproduced; a technique known as fast-mapping. By fitting these higher dimensional solution ensembles with polynomials, the original fast-mapping technique is extended to include T1 at arbitrary times. Accuracy of the 7th order polynomial is such that a minimum 96 per cent of the models are within a ±3 per cent relative error. This permits the rapid reproduction of full-Bloch solutions with a matrix multiplication and opens up surface NMR to time-series based inversion of single and multipulse data.
Wenke Sun,
Published: 11 August 2021
Geophysical Journal International, Volume 227, pp 2079-2094; https://doi.org/10.1093/gji/ggab320

Abstract:
SUMMARY Sources of expansion, such as nuclear explosions and magma chambers, can cause various surface effects, such as displacement and gravity field change on the Earth. The deformation simulation of an expansion source can be traced back to certain classic works nearly 70 yr ago, which were primarily based on a half-space earth model. However, the deformation of expansion sources in a spherical earth model has not been presented in the context of the semi-analytical viscoelastic Love numbers framework. In this study, we present the deformation Green's function for an expansion source in a spherical, non-rotating, viscoelastic, isotropic layered earth model using a semi-analytical method. We further obtain the analytical Green's functions for an elastic model in the near field. In addition, the influence of the stratification, curvature, and viscosity on the deformation is analysed using a layered spherical earth model. The suitability and effectiveness of the formulas are verified by numerical experiments and compared with previous results. The method presented in this paper is useful for explaining the viscoelastic deformations caused by the deep expansion source, for example, large magma chambers and upwelling of a mantle plume, in the layered viscoelastic earth model.
, Zachary E Ross, Kamyar Azizzadenesheli, Jack B Muir
Published: 11 August 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab309

Abstract:
We introduce a scheme for probabilistic hypocenter inversion with Stein variational inference. Our approach uses a differentiable forward model in the form of a physics informed neural network, which we train to solve the Eikonal equation. This allows for rapid approximation of the posterior by iteratively optimizing a collection of particles against a kernelized Stein discrepancy. We show that the method is well-equipped to handle highly multimodal posterior distributions, which are common in hypocentral inverse problems. A suite of experiments is performed to examine the influence of the various hyperparameters. Once trained, the method is valid for any seismic network geometry within the study area without the need to build travel time tables. We show that the computational demands scale efficiently with the number of differential times, making it ideal for large-N sensing technologies like Distributed Acoustic Sensing. The techniques outlined in this manuscript have considerable implications beyond just ray-tracing procedures, with the work flow applicable to other fields with computationally expensive inversion procedures such as full waveform inversion.
Published: 10 August 2021
Geophysical Journal International; https://doi.org/10.1093/gji/ggab318

Abstract:
Summary Gravity gradient tensors (GGTs) are used to investigate the density of subsurface structures in the Earth's crust, and can reduce ambiguity during data interpretation. However, costs and research area restrictions often prevent their application during surveys, thereby limiting their utility. To address this limitation, a matrix equation based on the Taylor series expansion that uses the gravity vector and its neighbors was formulated to obtain the GGTs. Higher-order derivatives of the gravity vector were utilized to constrain the calculation, which improved the accuracy of the transformation. Synthetic data were used to demonstrate that the proposed approach improved accuracy when the radial component of the gravity vector was transformed into GGTs. This approach was also applied to gravity data from the Otway Basin in Australia. Compared with the measured GGT, the results obtained using the proposed approach had a relative error of 0.46.
Published: 9 August 2021
Geophysical Journal International, Volume 227, pp 1893-1904; https://doi.org/10.1093/gji/ggab317

Abstract:
Summary A general-purpose model combining concepts from rational continuum mechanics, fracture and damage mechanics, plasticity, and poromechanics is devised in Eulerian coordinates, involving objective time derivatives. The model complies with mass, momentum, and energy conservation as well as entropy inequality and objectivity. It is devised to cover many diverse phenomena, specifically rupture of existing lithospheric faults, tectonic earthquakes, generation and propagation of seismic waves, birth of new tectonic faults, or volcanic activity, aseismic creep, folding of rocks, aging of rocks, long-distance saturated water transport and flow in poroelastic rocks, melting of rocks and formation of magma chambers, or solidification of magma.
Published: 9 August 2021
Geophysical Journal International, Volume 227, pp 2016-2043; https://doi.org/10.1093/gji/ggab310

Abstract:
SUMMARY In recent years, many higher-order and optimized schemes have been developed to reduce the dispersion error with the use of large grid spacing in finite-difference seismic waveform simulations. However, there are two problems in the application of coarse grids for heterogeneous media: the generation of artefact diffraction from the stair-step representation of non-planar interfaces and the inaccuracy of the calculated waveforms due to the interface error. Several equivalent medium parametrization approaches have been proposed to reduce the interface error of the finite-difference method in heterogeneous media. However, these methods are specifically designed for the standard (2,4) staggered-grid scheme and may not be accurate enough for coarse grids when higher-order and optimized schemes are used. In this paper, we develop a tilted transversely isotropic equivalent medium parametrization method to suppress the interface error and the artefact diffraction caused by the staircase approximation under the application of coarse grids. We use four models to demonstrate the effectiveness of the proposed method, and analyse the accuracy of each seismic phase related to the interface. The results show that our method can be used with higher-order and optimized schemes at 3 points per wavelength and produce satisfactory results.
Gerrit Hein, , , , György Hetényi, Rafael Abreu, Ivo Allegretti, Maria-Theresia Apoloner, Coralie Aubert, Simon Besançon, et al.
Published: 9 August 2021
Geophysical Journal International, Volume 227, pp 1996-2015; https://doi.org/10.1093/gji/ggab305

Abstract:
SUMMARY To constrain seismic anisotropy under and around the Alps in Europe, we study SKS shear wave splitting from the region densely covered by the AlpArray seismic network. We apply a technique based on measuring the splitting intensity, constraining well both the fast orientation and the splitting delay. Four years of teleseismic earthquake data were processed, from 723 temporary and permanent broad-band stations of the AlpArray deployment including ocean-bottom seismometers, providing a spatial coverage that is unprecedented. The technique is applied automatically (without human intervention), and it thus provides a reproducible image of anisotropic structure in and around the Alpine region. As in earlier studies, we observe a coherent rotation of fast axes in the western part of the Alpine chain, and a region of homogeneous fast orientation in the Central Alps. The spatial variation of splitting delay times is particularly interesting though. On one hand, there is a clear positive correlation with Alpine topography, suggesting that part of the seismic anisotropy (deformation) is caused by the Alpine orogeny. On the other hand, anisotropic strength around the mountain chain shows a distinct contrast between the Western and Eastern Alps. This difference is best explained by the more active mantle flow around the Western Alps. The new observational constraints, especially the splitting delay, provide new information on Alpine geodynamics.
, A M Lontsi, K Kremer, P Bergamo, , , D Fäh
Published: 9 August 2021
Geophysical Journal International, Volume 227, pp 1857-1878; https://doi.org/10.1093/gji/ggab314

Abstract:
Summary Single-station and array ambient vibration techniques are widely used in onshore environments, in particular to retrieve the subsurface structure and shear-wave velocity profiles. We apply these techniques offshore in Lake Lucerne (Switzerland) using single-station and array Ocean Bottom Seismometer (OBS) data. This lake has experienced tsunamigenic subaquatic slope failures in the past and still has sediment-charged slopes that might fail in the presence of a seismic or aseismic trigger. The application of traditional onshore methods offshore brings additional challenges related to the processing of recorded data. To overcome these challenges, we perform multibeam bathymetry surveys to precisely locate the OBS on the lake floor and airgun shootings to determine the orientation of the horizontal components of the seismometer and to correct the time drift of the recorder. Then we obtain surface-wave phase velocity dispersion curves of Scholte and Love waves, and Scholte wave ellipticity curves at six subaquatic slopes. After the estimation of the dispersion curves, we deal with their modal identification using mode attribution analysis. The shear-wave velocity and thickness of the sedimentary layers at the investigated slopes are inferred using a transdimensional Bayesian inversion algorithm. The resolved velocity profiles show very low shear-wave velocities in shallow lake sediments and allow us to improve the understanding of the local stratigraphy. This research contributes to the assessment of stability and tsunamigenic potential of subaquatic slopes in Lake Lucerne.
, Z Eilon, , P Audet
Published: 9 August 2021
Geophysical Journal International, Volume 227, pp 1879-1892; https://doi.org/10.1093/gji/ggab315

Abstract:
Summary Measurements of various physical properties of oceanic sediment and crustal structures provide insight into a number of geologic and geophysical processes. In particular, knowledge of the shear-wave velocity (VS) structure of marine sediments and oceanic crust has wide ranging implications from geotechnical engineering projects to seismic mantle tomography studies. In this study we propose a novel approach to non-linearly invert compliance signals recorded by co-located ocean-bottom seismometers and high-sample-rate pressure gauges for shallow oceanic shear-wave velocity structure. The inversion method is based on a type of machine learning neural network known as a mixture density neural network. We demonstrate the effectiveness of the mixture density neural network method on synthetic models with a fixed deployment depth of 2015 m and show that among 30,000 test models, the inverted shear-wave velocity profiles achieve an average error of 0.025 km/s. We then apply the method to observed data recorded by a broadband ocean-bottom station in the Lau basin, for which a VS profile was estimated using Monte Carlo sampling methods. Using the mixture density network approach, we validate the method by showing that our VS profile is in excellent agreement with the previous result. Finally, we argue that the mixture density network approach to compliance inversion is advantageous over other compliance inversion methods because it is faster and allows for standardized measurements.
, Yaoguo Li
Published: 9 August 2021
Geophysical Journal International, Volume 227, pp 2058-2078; https://doi.org/10.1093/gji/ggab316

Abstract:
SUMMARY Effective quantitative methods for integrating multiple inverted physical property models are necessary to increase the value of information and advance interpretation further to produce interpretable geology models through geology differentiation. Geology differentiation is challenging in greenfield exploration areas where specific a priori geological information is scarce. The main problem is to identify geological units quantitatively with appropriate 3-D integration of these models. The integration of multiple sources of information has been conducted with different unsupervised machine learning methods (e.g. clustering), which can identify relationships in the data in the absence of training information. For this reason, we investigate the performance of five different clustering methods on the identification of the geological units using inverted susceptibility, density, and conductivity models that image a synthetic geological model. We show that the correlation-based clustering yields the best results for the geology differentiation among those investigated by identifying the correlation between physical properties diagnostic of each unit. The result of the differentiation is a quasi-geology model, which is a model that represents the geology with inferred geological units and their spatial distribution. The resulting integrated quasi-geology model demonstrates that individually inverted models with minimal constraints have sufficient information to jointly identify different geological units.
Rongxin Fang, Jiawei Zheng, , Huanghui Lv, Chuang Shi, Jingnan Liu
Published: 6 August 2021
Geophysical Journal International, Volume 227, pp 1846-1856; https://doi.org/10.1093/gji/ggab313

Abstract:
Summary High-rate Global Navigation Satellite System (GNSS) has emerged as an effective method to recover seismic waveforms without saturation and drifts, but it has the limitation of relatively lower sampling rate and higher noise level compared to seismic instruments. In this study, we present a new seismogeodetic method by integrating GNSS and accelerometer data to obtain optimal real-time seismic waveforms. Unlike traditional integration methods based on GNSS techniques of relative positioning or precise point positioning, the new method uses a GNSS time difference technique and inherits its unique advantage in real-time and high-accuracy velocity solutions. Furthermore, by incorporating the tightly coupled structure, it can overcome the cascading problem and provide more accurate and robust waveforms compared to its loosely coupled counterpart. The performance of this method is first compared with the traditional loosely coupled approach in challenging environments through a set of shake table experiments. With three GNSS satellites, this approach method can improve the accuracy of velocities and displacements by 42 per cent and 87 per cent, respectively. With four or more GNSS satellites, the average improvements of the method reach 25 per cent and 41 per cent for the velocities and displacements, respectively. We then validate the full performances of the method through simulated shake table experiments and collocated GNSS and accelerometer data during the 2016 Mw6.6 central Italy earthquake. The simulated and real-event analyses demonstrate that the new integration method can take full advantage of the complementary characteristics of GNSS and accelerometer sensors. By providing more accurate and broadband velocity and displacement waveforms in a real-time or near-real-time manner, this method is quite promising in earthquake early warning and rapid source inversion.
Page of 343
Articles per Page
by
Show export options
Select all