Geophysical Journal International

SummaryThis article presents a new methodology to assist geophysicists in determining the first break event in a 3D seismic dataset using the well-known technique called Dynamic Time Warping algorithm (DTW), which is usually used to find the optimal alignment between two-time series. We used the optimal path from the cost matrix to identify the first-break in the seismogram using a few picks (seeds) made by an interpreter as a reference to perform this task. Furthermore, the data was pre-conditioned by the topographic and linear moveout to improve the method’s accuracy. To demonstrate the technique’s robustness, first, we applied the methodology in a synthetic seismic data. After demonstrating the efficiency of the algorithm, we applied the aforementioned methodology in the Polo-Miranga 3D seismic cube located in the Recóncavo sedimentary basin, Bahia-Brazil, and in the seismic data acquired from the Blackfoot field in Alberta, Canada. The high-quality results showed consistency in determining the first-break in all ranges of offsets, demonstrating an alternative way to accelerate this seismic processing step. Furthermore, we compared the results obtained by the proposed methodology with an algorithm based on comparing the short-time averages (STA) with long-time averages (LTA). Finally, we performed the static correction calculation to ensure that the time distortion resulting from the terrain and the low-velocity layer (LVL) was mitigated in shoot gathers and in the stacked section.

SummaryThis study presents the findings of a splitting analysis conducted on core-refracted teleseismic shear waves (SKS, SKKS, and PKS, called together as XKS) and local shear waves, obtained from a dense seismological network spanning the Kamchatka Peninsula. The objective of the study is to examine the pattern of mantle flow beneath the study area through the investigation of seismic anisotropy. The peninsula is situated at the northeastern end of the Kuril-Kamchatka subduction zone, where the Kuril trench intersects with the western boundary of the Aleutian trench. The dataset utilized in this study comprises waveform data from a dense network of seismic stations (99 broadband and short period stations for the local shear wave splitting analysis and 69 broadband stations for the SKS splitting analysis). The seismograms were downloaded from publicly available data repositories including the IRIS Data Management Center and the GFZ Data Services (GEOFON program). The dense station coverage allows us to investigate the lateral variations in anisotropy, providing insights into the flow patterns within the mantle. The processing of the combined data sets of local shear wave and teleseismic XKS waves allowed us to partially decipher the source of anisotropy in the mantle. Small delay (splitting) times (∼0.35 sec) observed from the local-S data suggest that anisotropy in the mantle wedge is relatively weak with lateral variations. Larger splitting times (∼1.1 sec) observed for the XKS waves relative to local S suggest that the main part of splitting on the XKS waves occurs in the subslab mantle. On the other hand, the rotational pattern of seismic anisotropy observed by both the local S and XKS waves suggests the presence of a toroidal flow at the NE edge of the subducting slab, which affects both the mantle wedge and subslab mantle. For the regions away from the edge of the slab, the mantle flow seems to be governed mainly by the drag of the lithospheric plate over the underlying asthenosphere.

SummaryPlate-coupling estimates and previous seismicity indicate that portions of the Makran megathrust of southern Pakistan and Iran are partially coupled and have the potential to produce future magnitude 7 + earthquakes. However, the GPS observations needed to constrain coupling models are sparse and lead to an incomplete understanding of regional earthquake and tsunami hazard. In this study, we assess GPS velocities for plate-coupling of the Makran subduction zone with specific attention to model resolution and the accretionary prism rheology. We use finite element model-derived Green's functions to invert for the interseismic slip deficit under both elastic and viscoelastic Earth assumptions. We use the model resolution matrix to characterize plate-coupling scenarios that are consistent with the limited spatial resolution afforded by GPS observations. We then forward model the corresponding tsunami responses at major coastal cities within the western Indian Ocean basin. Our plate-coupling results show potential segmentation of the megathrust with varying coupling from west to east, but do not rule out a scenario where the entire length of the megathrust could rupture in a single earthquake. The full subduction zone rupture scenarios suggest the Makran may be able to produce earthquakes up to Mw 9.2. The corresponding tsunami model from the largest earthquake event (Mw 9.2) estimates maximum wave heights reaching 2 to 5 meters at major port cities in the Northern Arabian Sea region. Cities on the west coast of India are less affected (1-2 m). Coastlines bounding eastern Africa, and the Strait of Hormuz, are the least affected (< 1 m).

SummaryUsing a Bayesian approach we compare anecdotal tsunami runup observations from the 29 December 1820 Flores Sea earthquake with close to 200,000 tsunami simulations to determine the most probable earthquake parameters causing the tsunami. Using a dual hypothesis of the source earthquake either originating from the Flores Thrust or the Walanae/Selayar Fault, we found that neither source perfectly matches the observational data, particularly while satisfying seismic constraints of the region. Instead both posteriors have shifted to the edge of the prior indicating that the actual earthquake may have run along both faults.

AbstractWe introduce a novel probabilistic framework for the solution of non-linear geophysical inverse problems in complex variables. By using complex probability distributions, this approach can simultaneously account for individual errors of real and imaginary data parts, independently regularize real and imaginary parts of the complex model, and still take into account cross-sensitivities resulting from a complex forward calculation. The inverse problem is solved by means of optimization. An application of the framework to complex resistivity (CR) imaging demonstrates its advantages over the established inversion approach for CR measurements. We show that CR data, with real and imaginary parts being subject to different errors, can be fitted adequately, accounting for the individual errors and applying independent regularization to the real and imaginary part of the subsurface conductivity. The probabilistic framework itself serves as a basis for the future application of global sampling approaches, such as Markov chain Monte Carlo methods.

SummaryThis study acquires the co-seismic deformation field of the MW7.4 earthquake that occurred in Maduo, China, on May 22, 2021, based on the BeiDou Navigation Satellite System (BDS), and the comparison with the results obtained by the Global Positioning System (GPS) reveals that the two systems are certain differences in their ability to acquire the horizontal co-seismic deformation field. The maximum difference in the horizontal co-seismic deformation is < 5 mm, and the maximum difference in the vertical co-seismic deformation is 8.7 mm. The dynamic displacement waveforms of the 2021 MW 7.4 Maduo earthquake acquired by BDS and GPS are very similar, which confirms that BDS can acquire ground-shaking images with an accuracy comparable to that of GPS. Based on the empirical relationship equation of the peak ground displacement (PGD) and moment magnitude (MW), this study verifies and calculates both the MW of the 2021 MW 7.4 Maduo earthquake and the error and finds that the MW can be quickly and accurately obtained by using the empirical PGD and MW equations, and this MW value can be used as a supplementary means of calibrating the MW of the large earthquake early warning (EEW) systems, which can be quickly determined by seismic wave data. Finally, by comparing the slip distributions inverted from the BDS and GPS co-seismic deformation fields, this study finds that BDS is equally effective as GPS.

SummaryHyuga-nada, off the Pacific coast of Kyushu along the Nankai Trough in southwest Japan, is one of the most active slow earthquake regions around Japan. We estimated the energies of shallow tremors and moments of shallow very low frequency earthquakes (VLFEs) in Hyuga-nada using data from a permanent onshore broadband network and temporary ocean bottom seismometer observations. The energies and moments of these slow earthquakes have a similar along-strike variation and are generally higher south of the subducted Kyushu-Palau Ridge than near the top of the ridge. This spatial variation is also related to the characteristics of slow earthquake migration. The along-strike migration speed was faster at initiation in the south, where the moments of slow earthquakes are higher. After migration entered the subducted Kyushu-Palau Ridge, its speed was decelerated with a parabolic pattern and their moments became smaller. Assuming a constant patch size of slow earthquakes, we estimated that the stress drop of VLFEs in the south of the subducted ridge was approximately three times higher than that near the top of the subducted ridge. According to our observations and a physical model, this stress drop difference between adjacent regions may cause parabolic migration. We also estimated the scaled energy of slow earthquakes from the ratio of the seismic energy rates of tremors to the seismic moment rates of accompanying VLFEs. The spatial variation in scaled energy is not identified inside the Hyuga-nada. Since the range of scaled energy is similar between the south and near the top of the subducted ridge, the apparent stress may be similar if the rigidity is the same. The dominant range of scaled energy of slow earthquakes in Hyuga-nada is 10−11.5–10−8.5. In addition to having similar or one order smaller values compared to other slow earthquake regions, the range of scaled energy in Hyuga-nada is broader. This broader range suggests wide range of characteristic time and various spectral features of slow earthquakes in Hyuga-nada. Based on a Brownian slow earthquake model, the wide range of characteristic time in this area suggests width variations of slow earthquake source area.

SummaryIn seismic tomography practices the Earth's surface is sometimes assumed to be either spherical or flat for convenience in forward modeling calculations. The effect of irregular surface topography on seismic wave propagation is thus ignored, resulting in biases in the phases and amplitudes of synthetic seismograms, which contribute to the residuals that are mapped into velocity structures in tomography inversions. In this study, we conduct a series of inversion experiments based on the adjoint waveform tomography method to quantitatively assess the topography effect on waveform-based inversion results. We first employ models with simplified topography to better highlight and quantify the topography effect. Results show that when topography effect is ignored in the forward modeling, it is mapped into velocity perturbations, leading to spurious velocity anomalies in tomography models. The strength of the spurious velocity anomaly is quasi-linearly related to locally averaged topography gradient. Our inversion experiments demonstrate that in places of strong topography variation, such as the Longmenshan Fault Zone region where the 2008 Mw7.9 Wenchuan earthquake occurred, topography effect can lead to spurious relative velocity anomalies of up to 10 per cent, which cannot be ignored in waveform-based tomography inversions.

SummarySeismic anisotropy is key to constrain mantle flow, but it is challenging to image and interpret it. Existing large-scale tomography models of seismic anisotropy typically show large discrepancies, which can lead to completely distinct geodynamical interpretations. To better quantify the robustness of anisotropy tomography, we create a 2-D ridge-to-slab geodynamic model and compute the associated fabrics. Using the resulting 21 elastic constants we compute seismic full waveforms, which are inverted for isotropic and radially anisotropic structure. We test the effects of different data coverage and levels of regularisation on the resulting images and on their geodynamical interpretation. Within the context of our specific imposed conditions and source-receiver configuration, the retrieved isotropic images exhibit substantial artificial slab thickening and loss of the slab’s high velocity signature below ∼100 km depth. Our results also show that the first order features of radial anisotropy are well retrieved despite strong azimuthal anisotropy (up to 2.7 %) in the input model. On the other hand, regularisation and data coverage strongly control the detailed characteristics of the retrieved anisotropy, notably the depth-age dependency of anisotropy, leading to an artificial flat depth-age trend shown in existing anisotropy tomography models. Greater data coverage and additional complementary data types are needed to improve the resolution of (an)isotropic tomography models.

SummaryIn the context of Bayesian inversion, we consider sequential Monte Carlo (SMC) methods that provide an approximation of the posterior probability density function and the evidence (marginal likelihood). These particle approaches build a sequence of importance sampling steps between gradually-tempered distributions evolving from the prior to the posterior PDF. To automate the definition of the tempering schedule, adaptive sequential Monte Carlo (ASMC) allows tuning the temperature increments on-the-go. One general challenge in Bayesian inversions is the computational burden associated with expensive, high-fidelity forward solvers. Lower-fidelity surrogate models are interesting in this context as they can emulate the response of expensive forward solvers at a fraction of their cost. We consider surrogate modeling within ASMC and introduce first an approach involving surrogate modeling only, in which either prior samples are used to train the surrogate, or the surrogate model is re-trained by updating the training set during the inversion. In our implementation, we rely on polynomial chaos expansions for surrogate modeling, principal component analysis for model parametrization and a ground-penetrating radar cross-hole tomography problem with either an eikonal or finite-difference time-domain solver as high-fidelity solver. We find that the method based on re-training the surrogate during the inversion outperforms the results obtained when only considering prior samples. We then introduce a computationally more expensive multifidelity approach including a transition to the high-fidelity forward solver at the end of the surrogate-based ASMC run leading to even more accurate results. Both methods result in speed-ups that are larger than one order of magnitude compared to standard high-fidelity ASMC inversion.

SummaryDense array observation and seismic interferometry have revolutionized the imaging schemes of the earth structure. It is becoming possible to directly obtain the lateral variation of the earth's structure by applying array-based methods such as the cross-correlation beamforming (CBF) of the ambient noise to the subsets of the dense array, without tomography. CBF has been proven to extract the azimuth-averaged multi-mode surface wave dispersion curves. However, the resolution of the dispersion image generated by conventional CBF is low at high frequencies in the frequency-velocity (f-v) domain. Moreover, the irregular array geometry and uneven source distribution would bias the result of CBF, especially for the estimation of azimuth-dependent velocity. In this paper, two beamforming schemes are suggested to improve the resolution of multi-mode dispersion images in the f-v domain. Firstly, the geometrical spreading of the wavefield is corrected to enhance the amplitude at high frequency (or large distance) and thereby improve the resolution of the dispersion image at high frequency. We call this scheme weighted correlation beamforming (WCBF). The azimuth-averaged velocity can be estimated with sufficient resolution using WCBF by stacking the beamforming output at each azimuth. We show that WCBF is the 2D Fourier transform of the spatial wavefield from the viewpoint of the wavefield transform. Secondly, a modified beamforming scheme (MCBF) is suggested to reduce the effect of uneven source and/or irregular array geometry. The delay and summation in MCBF are performed only for plane waves incident from the stationary phase region. The azimuth-dependent velocity can therefore be estimated by MCBF with less dependence on the array geometry, as well as on the uneven source distribution. In terms of the estimation of azimuth-averaged phase velocity, we show the F-J method, another array-based method for extracting multi-mode surface waves from ambient noise using the Fourier-Bessel transform, is the azimuth-averaged version of WCBF. The reliability of WCBF and MCBF is verified based on the synthetic and field data using the array with different geometry. The dispersion image of multi-mode Rayleigh wave phase velocity at local and regional scales can be generated by WCBF or MCBF with high resolution. In particular, multi-mode dispersion curves at the local scale can be measured by MCBF with sufficient accuracy using quite short recordings from hours to days. This offers the possibility of a rapid assessment of the media properties.

SummaryTo understand the seismic hazard of a subduction zone, it is necessary to know the geometry, location, and mechanical characteristics of the interplate boundary below which an oceanic plate is thrust downward. By considering the azimuthal dependence of converted P-to-S (Ps) amplitudes in receiver functions (RFs), we have detected the interplate boundary in the Makran subduction zone, revealing significant seismic anisotropy at the base of the accretionary wedge above the slab before it bends down beneath the Jaz Murian basin. This anisotropic feature aligns with a zone of reduced seismic velocity and a high primary/secondary wave velocity ratio (Vp/Vs), as documented in previous studies. The presence of this low-velocity highly anisotropic layer at the base of the accretionary wedge, likely representing a low-strength shear zone, could possibly explain the unusually wide accretionary wedge in Makran. Additionally, it may impact the location and width of the locked zone along the interplate boundary.

SummaryA major earthquake shook the Chinese county of Maduo, located in the Songpan-Ganzi terrane on the Tibetan Plateau, on 21 May 2021. Here, we investigate the postseismic deformation process of this event, with the aim to understand the fault geometry, friction behaviour and regional rheology. To keep the self-consistency between co- and post-seismic deformation models, we first constrain the fault geometry and coseismic slip model of this event, which are directly used in modelling the postseismic deformation. The coseimsic slip model reveals that the majority of coseismic slip is confined at the middle (3-15 km) of the brittle layer, leading to significant shallow slip deficit. Secondly, we obtain the postseismic deformation in the first 450 days following the 2021 Maduo earthquake using the GPS and InSAR displacement time series data. Thirdly, a combined model incorporating afterslip and viscoelastic relaxation is built to explain the observed postseismic deformation. Our results suggest that the viscoelastic relaxation effect should be considered in the observation period, in order to avoid the unphysical deep afterslip in the ductile lower crustal layer. Combined analysis on viscosities inferred from this study and previous studies suggests a weak lower crust with steady-state viscosity of 1018 ∼ 1019 Pa s beneath the Songpan-Ganzi terrane, which may give rise to the distributed shear deformation and the development of subparallel secondary faults within the terrane. Besides, the inferred afterslip on uppermost patches of the middle fault segment suggests a rate-strengthening frictional behaviour that may be related to the coseismic slip deficit and rupture arrest of the Maduo earthquake.

SummaryAmbient noise technology can efficiently extract surface wave signals from seismic background noise and has been extensively utilized in imaging lithospheric structures. However, retrieving crustal body wave signals, such as PmP or SmS phases, still poses a challenge. Only a limited number of reports have successfully extracted these regional-scale body wave signals from ambient noise in only a few limited study areas. It remains unclear why these signals are difficult to retrieve from ambient noise data.To investigate the mechanism of recovering body wave signals in noise cross-correlations, we calculate cross-correlation functions at four regions and observe the similarity of the recovered body waves. Through a series of synthetic simulations, we demonstrate that the appearance of body wave signals in noise cross-correlations is closely related to the distribution of noise sources. Among these signals, the post-critical SmS wave proves to be the most readily recoverable from ambient noise data, primarily stemming from distant sources. In contrast, the recovery of P-wave requires the array to be in proximity to the sources. Our experiments also reveal that the main origin of PL waves is the multiple reflections of S-waves propagating in the crust.

SummaryFor a more quantitative discussion of slow earthquake activity, we evaluated the detectable limits of very low frequency earthquakes (VLFEs), which are seismic slow earthquakes observed in very low frequency (< 0.05 Hz) bands in the Nankai subduction zone. We performed numerical simulations using a local three-dimensional model and used the observed noise level of permanent broadband seismometers. First, we investigated the effects of the source-time functions on the maximum amplitudes of the VLFE signals at a certain station. The maximum amplitudes of the VLFE signals were controlled by the VLFE moment rate. The detectable limit of VLFEs at each source location can be defined as the lowest moment rate of detectable VLFEs, which radiate signals larger than the noise levels of any component at ≥ 3 stations. For inland seismometers only, the detectable limits of VLFEs at deep (30–40 km) and shallow (≤ 10 km) depths were 1012–1012.3 and 1012.7 Nm/s, respectively. Due to the geometrical spreading of VLFE signals and large noise levels in horizontal components, offshore seismometers improved the detectability of shallow VLFEs in regions where seismometers were densely deployed. Based on our detectability and published catalogs, shallow slow earthquakes are less active south-southwest off the Kii Peninsula, where geodetic studies expect mechanical coupling.

SummaryWe have derived an analytical approximate expression to estimate the delay in diving seismic waves due to thin layers of CO2. The expression is valid for high frequencies and can be used to estimate the delay in diving waves at seismic frequencies for large separations between the source and receiver (offset). The approximation may be used to assess CO2 detection limits using diving waves and to support survey planning for CO2 monitoring and full-waveform inversion (FWI) cycle skipping analysis. In this study, we analyze the diving-wave response to a thin layer of CO2 for band-limited data using acoustic finite-difference modeling, and compare the results against the analytical calculations. We find that the responses are offset-dependent and related to double- and single-leg interactions between the diving waves and the CO2. To test the methods, we created a synthetic representation of the 2010 subsurface conditions for the top CO2 layer at the Sleipner storage complex in the North Sea, by combining base and monitor post-stack seismic data with field velocity trends. Using the acoustic finite-difference method, we model pre-stack data that captures the complexity of field data and demonstrate the use of the diving-wave delay for CO2 migration monitoring and CO2 thin layer detection.

SummaryMagnetotelluric data from Mount Tongariro has been analysed using an unstructured tetrahedral finite-element inversion code that incorporates topography, which was not included in previous analysis of these data. Incorporating topography adds information, which stabilises the resistivity inversion modelling, and for the first time allows details of the shallow hydrothermal system and its relationship with the underlying magmatic system to be resolved. Specifically, an electrically conductive zone between 4 and 12.5 km depth marks the underlying magmatic system, which is shown to directly connect via conductive pathways to the area where the most recent phreatic eruptions at Tongariro occurred in 2012. The resultant phreatic eruptions in August and November 2012 showed no new magmatic component to the eruption deposits. Nevertheless, by combining the magnetotelluric resistivity image with relocated seismicity, we can see that seismicity (a proxy for magma ascent) migrated from the top of the magmatic system into the hydrothermal system in the months preceding these eruptions. Magmatic interaction with the extant hydrothermal system likely caused the over pressurisation for the phreatic eruption. This work highlights the utility of combining geophysical methods and petrological data to constrain phreatic eruption processes.

SummaryHydraulic Fracturing (HF) often stimulates the local earthquake productivity which provides a unique opportunity to characterize the crustal heterogeneities, reservoir properties, and fluid injection effects. However, the velocity models acquired solely based on the arrival time records are often undermined due to the seismic network coverage and interpolation techniques. Instead, we adopt the waveform-based approach to apprehend; (1) structural heterogeneities, (2) reservoir distribution, and (3) signatures of the injected fluid in the Weiyuan shale gas field (WSGF). We categorize the waveforms into dominant high and low frequencies based on the qualitative inspection and frequency index analysis of the seismic waveforms. We first inspect the waveform to access the potential controlling mechanisms (source, site, and path effects) at both single and multiple stations in different azimuthal orientations. As a result, we find the path effect as a dominant factor to influence the waveform characteristics, e.g., S-wave amplitude, and frequency. Subsequently, to localize the path effect, we conduct an in-depth examination of events within 10 km of each seismic station and classify the waveform records using their frequency indices. Notably, certain stations record a significant proportion of low-frequency waveforms (LFWs) (up to 20%), while others have limited occurrences (∼1%) indicating suspected anomalous zones. Afterward, we identify two suspected anomalous zones based on LFWs intensity and ray tracing map. Both zones are in close proximity to fault zones and preserved reservoirs with no HF activities, where fault damage zones or the fluid-rich reservoir may contribute to our observed LFWs.

SummaryThe south-eastern offshore of the Boso peninsula in Japan periodically experiences short-term slow slip events (SSEs) every few years. On June 2018, a SSE occurred with the maximum surface horizontal displacement reaching up to 4.7 cm by according to the operational global navigation satellite system (GNSS) network. This study performed a time series analysis of interferometric synthetic aperture radar (InSAR) with Sentinel-1 SAR images to investigate detailed spatial pattern of surface displacements caused by the SSE. With the assistance of an atmospheric delay correction with a regional numerical weather model output, the InSAR time series analysis successfully captured displacement signals in three paths, whose maximum amplitudes in line-of-sight directions were 1.46 cm, 1.86 cm and -0.80 cm. A checkerboard test revealed that the resolution of the slip inversion was higher when InSAR was used than that using GNSS, especially in and around the inland. The slip inversion with the actual displacement data derived from the InSAR time series analysis was performed with the L-curve optimization, showing that the estimated slip area was concentrated offshore south-eastward from the Boso peninsula with the maximum slip of 5 cm and the estimated moment magnitude of 6.4. As similar to previous SSEs in the Boso peninsula, a seismic swarm simultaneously occurred in the down-dip area adjacent to the estimated slip with the SSE occurrence, suggesting a different friction characteristics between them. This study demonstrates usefulnesses of the InSAR observation for capturing detailed spatial characteristics of small-displacement events like SSEs and of the hybrid use of the externally derived delay correction with the time series analysis to improve the displacement detection accuracy.

SummaryLow-frequency earthquakes (LFEs) are seismic phenomena with the shortest timescale among various slow earthquakes observed on broadband time scales. To understand the nature of such a broadband slow phenomenon, it is important to investigate the rupture evolution process of individual slow events, such as LFEs. Here, we investigated the moment-duration relationship of LFEs at plate interfaces and volcanic regions, and showed that the moment-duration relationship of both tectonic and volcanic LFEs is characterised by a moment proportional to the cubic duration, similar to that in ordinary earthquakes. The difference between our obtained moment-duration relationship and the broadband scaling suggests that the evolution process of LFEs may not be controlled, but only triggered by the slow earthquakes with longer durations, such as slow slip events driven by aseismic diffusion. The seismic moments of the LFEs are approximately three orders of magnitude smaller than those of ordinary earthquakes with similar durations. This result indicates that LFEs have rupture growth similar to that of ordinary earthquakes, although the rupture velocity and/or stress drop are much smaller. Considering the hypocentre spread of LFEs, the estimated rupture velocity and stress drop were approximately 100 m/s–1 km/s and 2 kPa–1 MPa, respectively. Additionally, the estimated moment magnitudes are much larger than the local magnitudes determined based on the maximum amplitudes, which is due to the longer durations and resultant smaller amplitudes of LFEs than those of ordinary earthquakes.