Geophysical Journal International

SummaryMagneto-Coriolis (MC) modes in Earth’s fluid core involve oscillations sustained by the combined effect of the Lorentz and Coriolis forces. Here, we investigate the properties of MC modes that involve purely axisymmetric flow, which we term axiMC modes. We provide a basic description of the wave dynamics of these modes, and simple predictions for the expected scalings of their frequency ω, decay rate λ, and quality factor Q based on a uniform ambient magnetic field. In particular, Q scales with the Elsasser number Λ, which depends on the square of the r.m.s. strength of the azimuthally averaged meridional field. When Λ > 1, Q > 1 and axiMC modes may be excited; when Λ ≪ 1, Q ≪ 1 and axiMC modes revert to quasi-free magnetic decay modes. We present computations of axiMC modes in an inviscid, electrically conducting sphere for two idealized ambient magnetic field configurations, a uniform axial field and an axial poloidal field. We show that a flow gradient in the axial direction is a key property of axiMC modes. For the uniform axial field, ω, λ and Q follow the scalings expected for a uniform field. For the axial poloidal field, the structure of the modes changes substantially when Λ ≳ 1, becoming more concentrated in regions of lower field strength. The combination of this structural change and advection of field lines by flow significantly increases λ, resulting in a Q that remains close to 1 even at high Λ. For a magnetic field strength inside the Earth’s core of a few mT, the gravest axiMC modes are expected to have periods in the range of one thousand to a few thousand years and a Q not substantially above 1. AxiMC modes may be connected to a part of the observed millennial changes in Earth’s magnetic field, may exchange axial angular momentum with the mantle, and hence may also explain a part of the observed millennial changes in length of day.

SummaryMulti-channel analysis of surface wave (MASW) is a non-destructive technique to characterize the sub-surface using the dispersive nature of Rayleigh waves. Field dispersion curves are inverted to predict the shear wave velocity structure of the ground and pavement profile. Adjusting the dynamic properties of the initially assumed soil profile necessitates information regarding the dominant sensitive layers. Therefore, a swift and accurate computation of the Jacobian of phase velocity is essential to generate an appropriate shear wave velocity profile and accelerate the inversion process. This is especially crucial for the 2D MASW survey, which requires hundreds of 1D inversions to create a high resolution 2D profile. Available numerical methods are computationally expensive and often suffer from instabilities for highly sensitive layers. The existing analytical methods involve mathematical complexities and require rigorous treatment. Furthermore, they are time-consuming and often found to be marginally faster than the numerical methods. Based on the fast delta matrix algorithm, the paper presents a new efficient analytical formulation of the Jacobian matrix of modal phase velocities concerning the layer parameters. The proposed algorithm leverages the simpler and fewer matrix elements of the fast delta matrix, thus significantly reducing the number of mathematical operations required. Additionally, it reduces the algorithm's cost by factorizing non-zero elements, thereby markedly reducing the computational time. Five different types of synthetic earth models are adopted from the published literature to validate the accuracy and efficacy of the newly developed algorithm. The presented work will significantly benefit the practicing engineers and geophysicists in processing field MASW test data.

SummaryMonitoring volcanic deformation is crucial for understanding volcanic behavior, but challenges like limited GNSS coverage, infrequent SAR data acquisitions, and coherence loss during eruptions complicate this task. Our study on the 2018 eruption of Sierra Negra utilizes Sentinel-1A/B images to track surface deformation patterns that revert to their initial state across three phases: before, during, and after the eruption. We implemented an adaptive workflow using the shortest temporal baseline of consecutive SAR image pairs, including InSAR, optical flow, and pixel offset tracking methods, to accurately capture surface displacement linked to the dynamics of the magma reservoir in the caldera and a nearby (sub-) horizontal dike. Results show that while the caldera subsided gradually over two months during lava flow (initially at a rate of several meters) until it began to uplift again, the northern region alternated between uplift and subsidence twice in the line-of-sight (LOS) direction. This pattern suggests repeated magma injections into the sub-horizontal dike sustained the lava flows from the northeastern fissure. The one-day difference between SAR images from ascending and descending tracks enabled us to estimate the underground magma transfer rate at approximately 60 m/h, which aligns with the magma migration trajectory indicated by seismic data. By integrating InSAR and offset tracking methods, we provide a comprehensive view of surface displacement throughout the volcanic eruption cycle.

AbstractCrystallographic preferred orientation (CPO) of peridotite minerals is frequently invoked to explain the widespread dependence of seismic wave speed on propagation direction in Earth’s mantle — a property known as seismic anisotropy. As established by rock mechanics experiments, CPO constitutes a direct signature of past and ongoing strain regimes experienced by rocks during mantle flow. Therefore, an improved understanding of CPO generation promises to yield valuable information on the rheology and corresponding deformation mechanisms activated through mantle dynamics. Simulating CPO in geodynamical models is computationally challenging and has often been restricted to steady-state mantle flows. However, within Earth’s vigorously convecting mantle the steady-state assumption is questionable, thus motivating the need to couple CPO simulations with time-evolving mantle flow models. Here, we present a new Python implementation of the D-Rex CPO model, called PyDRex, which predicts salient features of mineral grain size and orientation evolution whilst providing a well-documented, user-friendly interface that supports flexible coupling to geodynamical modelling frameworks. PyDRex also packages numerous post-processing routines for strain analysis and visualisation of grain orientation distributions. We provide a set of benchmark simulations based on previous D-Rex implementations that validate PyDRex and demonstrate sensitivities to model parameters for both steady-state and time-dependent flows. Analysis of benchmark results highlights the role of dynamic recrystallisation in controlling competing grain growth in both the softest and hardest crystallographic orientations. When employing a commonly used value for the grain boundary mobility parameter (M* = 125), we also find that transient CPO textures are generally not well resolved if crystals are represented by fewer than 5000 ‘grains’ (weighted orientation samples) — a configuration rarely employed in most previously published studies. Furthermore, kinematic corner-flow models suggest that CPO produced at mid-ocean ridges has a non-linear dependence on depth, which implies that even ostensibly simple mantle flows can result in complex distributions of seismic anisotropy. Our analyses motivate further experimental calibration of parameters controlling dynamic recrystallisation and potential improvements to the numerical treatment of subgrain nucleation.

SummaryForward modeling is crucial for seismic data processing, which is the core of reverse time migration and full-waveform inversion. Numerical simulation based on conventional elastic wave equations in stationary solids neglects the fluidity of fluids (e.g., seawater), making it difficult to simulate the propagation of seismic waves in moving fluids accurately. To solve the problem, we start with classical equations of fluid mechanics and derive a new set of elastic wave equations that can be used to simultaneously model wave propagation both in moving fluids and stationary solids. For high-precision numerical simulations, a staggered-grid finite-difference scheme is used to solve the proposed equations. Numerical tests on a homogeneous uniformly moving model demonstrate that the dynamic and kinematic characteristics (e.g., wavelength, amplitude) of elastic waves in moving fluids are quite different from those in stationary medium. Forward modeling for a two-layer model that has a flowing water layer and a stationary rock layer is used to study the reflection and transmission patterns of elastic waves in the solid-fluid interface. With the help of the superposition principle of vectors and Snell's law, the transmission angles can be easily calculated. A further test for a more complex stratified model indicates that the energy and travel time differences of reflected waves are expected to be evidence for the identification of moving fluids. Numerical experiments on the Marmousi II model demonstrate that the relative wavefield error is positively correlated with the maximum moving velocity and the wavelet dominant frequency.

SummaryDispersion curves of surface waves are widely used for the inversion of subsurface structures. To extract dispersion curves, many methods have been developed. Among them, multichannel analysis of surface waves such as slant stack and frequency-Bessel transform can extract not only the fundamental mode but also overtones. Inversion with overtones is proven to be more stable and has better resolution at greater depths. However, with a limited number of array receivers, artifacts and misfits due to array-layout effects arise in the dispersion spectra and impede the identification of dispersion curves. We evaluate the array-layout effects in the frequency-Bessel transform and calculate the array response functions which can help to mitigate artifacts and calibrate dispersion curves. We apply this technique to synthetic simulated, active source and ambient noise data. The artifacts caused by array-layout effects can be mitigated, which helps the identification of dispersion curves. We further calculate the Pearson correlation coefficient between the array response function and the dispersion spectrum section. It is used to calibrate the biases produced by the array-layout effects if we select dispersion curves by maximum values. The confidence intervals of the dispersion curves are then determined based on the correlation coefficients. It is helpful for the design of array layouts according to the investigation depths of interest.

SummaryPowerful infrasound (acoustic waves <20 Hz) can be produced by explosive volcanic eruptions. The long-range propagation capability, over hundreds to thousands of kilometers, of atmospheric infrasound motivates the development of regional or even global scale volcano-infrasound monitoring systems. Infrasound propagation paths are subject to spatiotemporal atmospheric dynamics, which lead to deviations in the direction-of-arrival (back-azimuth) observed at sensor arrays and contribute to source location uncertainty. Here we further investigate the utility of empirical climatologies combined with 3-dimensional ray-tracing for providing first-order estimates of infrasound propagation paths and back-azimuth deviation corrections. The intended application is in scenarios requiring rapid or precomputed infrasound propagation calculations, such as for a volcano-infrasound monitoring system. Empirical climatologies are global observationally based function fitting models of the atmosphere, representing robust predictors of the bulk diurnal to seasonal atmospheric variability. Infrasound propagation characteristics have previously been shown to have strong seasonal and diurnal components. At the International Monitoring System (IMS) infrasound station IS22, New Caledonia, quasi-continuous multi-year infrasound array detections show oscillating azimuthal variations for arrivals from volcanoes in Vanuatu, including Yasur (∼400 km range), Ambrym (∼670 km range), and Lopevi (∼650 km range). We perform 3-dimensional ray-tracing to model infrasound propagation from the Ambrym and Yasur volcano locations to IS22 every six hours (00:00, 06:00, 12:00, and 18:00 UTC) for every day of 2004 and 2019 for Ambrym and Yasur, respectively and evaluate the results as compared to the multi-year observations. We assess a variety of models and parameterizations, including both empirical climatologies and hybrid descriptions; range-independent and range dependent atmospheric discretizations; and unperturbed and perturbed range-independent empirical climatologies. The hybrid atmospheric descriptions are composed of ERA 5 reanalysis descriptions from the European Centre for Medium-Range Weather Forecasts (ECMWF) below ∼80 km altitude combined with empirical climatologies above. We propose and employ simple parametric perturbations to the empirical climatologies, which are designed to enhance the stratospheric duct and compensate for missing gravity wave perturbations not included in the climatologies, and thereby better match observations. We build year-long back-azimuth deviation interpolations from the simulations and compare them with three different multi-year array detection datasets from IS22 covering from 2003 up to 2022. Through a systematic comparison, we find that the range-independent empirical climatologies can capture bulk azimuth deviation variability and could thus be useful for rapid infrasound propagation calculation scenarios, particularly during favorable sustained propagation ducting conditions. We show that the hybrid models better describe infrasound propagation during periods of weak stratospheric ducting and during transient atmospheric changes such as stratospheric wind reversals. Overall, our results support the notion that climatologies, if perturbed to compensate for missing gravity wave structure, can improve rapid low-latency and precomputed infrasound source discrimination and location procedures.

SummaryStress accumulation on the plate interface of subduction zones is a key parameter that controls the location, timing and rupture characteristics of earthquakes. The diversity of slip processes occurring in the megathrust indicates that stress is highly variable in space and time. Based on GNSS and InSAR data, we study the evolution of the interplate slip-rate along the Oaxaca subduction zone, Mexico, from October 2016 through October 2020, with particular emphasis on the pre-seismic, coseismic and post-seismic phases associated with the June 23, 2020 Mw 7.4 Huatulco earthquake (also known as La Crucecita earthquake), to understand how different slip regimes contribute to the stress accumulation in the region. Our results show that continuous changes in both the aseismic stress-releasing slip and the coupling produced a high stress concentration (i.e., Coulomb Failure Stress (CFS) of 80 kPa) prior to the event on the region with the highest moment release of the Huatulco earthquake (between 17 and 30 km depth) and a stress deficit zone in the adjacent updip region (i.e., shallower than 17 km depth with CFS around -90 kPa). This region under negative stress accumulation can be explained by possible recurrent shallow Slow Slip Events (SSE) offshore Huatulco as well as by the stress shadow from adjacent locked segments. Absent in the literature, the shallow rupture is characterized by a secondary slip patch (between 7 and 14 km depth) that overlaps with the highest concentration of aftershocks. Two months prior to the event, a Mw 6.6 long-term SSE also occurred about 80 km northwest from the hypocenter, between 25 and 55 km depth. Transient increments of the interplate coupling around the adjacent 1978 (Mw 7.8) Puerto Escondido rupture zone correlate with the occurrence of the last three SSEs in Oaxaca far downdip of this zone, possibly associated with along-dip fluid diffusion at the subduction interface. Throughout the four-year period analyzed, the interface region of the 1978 event experienced a high CFS build up of 80-150 kPa, primarily attributable to both the co-seismic and early post-seismic slip of the Huatulco rupture, that, considering the 55 year average return period of the region, indicates large earthquake potential near Puerto Escondido. Continuous monitoring of the interplate slip-rate thus provides a better estimation of the stress accumulation in seismogenic regions than those given by long-term, time-invariant coupling models, and improves our understanding of the megathrust mechanics where future earthquakes are likely to occur.

SummaryMonitoring networks increasingly aim to assimilate data from a large number of diverse sensors covering many sensing modalities. Bayesian optimal experimental design (OED) seeks to identify data, sensor configurations, or experiments which can optimally reduce uncertainty and hence increase the performance of a monitoring network. Information theory guides OED by formulating the choice of experiment or sensor placement as an optimization problem that maximizes the expected information gain (EIG) about quantities of interest given prior knowledge and models of expected observation data. Therefore, within the context of seismo-acoustic monitoring, we can use Bayesian OED to configure sensor networks by choosing sensor locations, types, and fidelity in order to improve our ability to identify and locate seismic sources. In this work, we develop the framework necessary to use Bayesian OED to optimize a sensor network’s ability to locate seismic events from arrival time data of detected seismic phases at the regional-scale. This framework requires five elements:(i) A likelihood function that describes the distribution of detection and travel time data from the sensor network,(ii) A prior distribution that describes a priori belief about seismic events,(iii) A Bayesian solver that uses a prior and likelihood to identify the posterior distribution of seismic events given the data,(iv) An algorithm to compute EIG about seismic events over a dataset of hypothetical prior events,(v) An optimizer that finds a sensor network which maximizes EIG.Once we have developed this framework, we explore many relevant questions to monitoring such as: how to trade off sensor fidelity and earth model uncertainty; how sensor types, number, and locations influence uncertainty; and how prior models and constraints influence sensor placement.

SummaryThe seismic wavefield is fully described by translation, rotation and strain. Until recently, the seismological community has not been able to measure rotation directly with portable sensors with satisfactory resolution. Portable blueSeis-3A (Exail) sensors allow measuring 3 components of rotational motions. Co-located with conventional seismometers, one can locally observe six degrees of freedom (6 DoF) of ground motion. To test the performance of the rotational sensors, an active source experiment was carried out in Fürstenfeldbruck, Germany, in November 2019. Five explosions with different yields and distances from our sensors, were fired. In a first stage, eight rotational sensors were deployed inside a bunker next to a seismometer. In a second stage, the rotational sensors were installed in two clusters of 4 sensors each. We compare the back azimuths derived using two different methods: (i) a method using horizontal rotational components and (ii) a standard polarization analysis using only 3C translational data. Back azimuths calculated using rotational data for 5 explosions have an average 10.2○ deviation from the theoretical back azimuths. Estimates using 3C translational data for the first stage of the experiment have a 1.8○ deviation from the theoretical back azimuth. For the second stage we found a 29.4○ deviation using the seismometer from stage 1. We conclude that within our distance range from 50 m to 1070 m, all rotational sensors provide reliably back azimuths of explosive sources when using only horizontal rotational components. For future applications of rotational sensors in other environments this is promising as back azimuths can be derived reliably.

SummaryA new global mascon solution using GRACE and GRACE Follow-On data is co-estimated with Satellite Laser Ranging (SLR) measurements to seven major geodetic satellites. This combined solution is compared with an otherwise similar GRACE-only solution to determine improvements in the estimate. We find similar performance between both solutions in the recovered mass change, but significant improvements in the associated errors in the combination solution. Errors in recovered basin mass change are 10-20% better for the combination solution across all basin sizes, with the greatest improvements in high latitude ice sheets. These results lead to our recommendation that all GRACE Level 3 mascon and spherical harmonic user-oriented gridded solutions should include SLR information during the solution inversion. As an ancillary contribution, we also provide validation of the choice of truncated spherical harmonics used in determining GRACE Technical Note 14, the current recommended mechanism for including SLR information with GRACE solutions in post-processing.

SummaryThe 1929 Grand Banks earthquake (Mw 7.1) generated a large landslide that broke submarine cables near the source area. A large tsunami was generated by the large landslide that killed 28 people on the southern coast of the Burin Peninsula, Canada. The tsunami was also observed at a tide gauge in Halifax, Canada. In this study, the initial landslide mass distribution was estimated by fitting the first pulse of the tsunami waveform observed at Halifax and the timings of the submarine cable break with those computed values. A deep-sea landslide tsunami computation method was developed by modifying a previously developed Tsunami Squares method with a different friction term and an effect of deep water. The results showed that the tsunami waveform and timing of the cable break were well explained by the tsunami computed from the initial landslide mass volume of about 450 km3 located along the continental slope near the source area. The method developed to compute deep-sea landslide tsunamis should be useful for studying other types of landslide tsunamis.

SummaryWe conduct a comprehensive comparison of Ice Mass Balance (IMB) estimates for Greenland derived from satellite observations of ice Surface Elevation Changes (SEC), gravity and GNSS observations. Our analysis integrates data from the ICESat and CryoSat-2 satellite altimetry missions, augmented by optical stereo-imagery for peripheral glaciers, and GRACE satellite gravimetry mission, spanning the 2003-2008 and 2011-2015 periods. We also consider three firn densification models (FDM) and five Glacial Isostatic Adjustment (GIA) models for correcting the datasets for these effects when necessary. Our results reveal significant differences among FDM corrections applied to SEC observations, with particularly large variations in IMB estimates reaching up to 90 Gt/yr. To address this, we develop an innovative method for estimating equivalent firn corrections to the ice elevation observations, based on a least-squares fit of filtered ice SEC observations to GRACE mass-change estimates. This approach is both simple and independent from climate models assumptions and shows minimal sensitivity to GIA model differences. Using this method, we estimate IMBs for Greenland at -217.6 ± 15.7 Gt/yr for 2003-2008 and -253.2 ± 18.8 Gt/yr for 2011-2015. Importantly, these values indicate an acceleration of the thinning rate, not consistently captured by the IMB estimates inferred from the ice SEC observations corrected by FDMs. Finally, we compute elastic ground deformation induced by ice mass change during 2011-2015, using the four proposed mass-variation distributions and compare the predicted vertical velocities with GNSS observations in Greenland, accounting for all GIA models. While all models are consistent with most of the GNSS-derived uplift rates, they cannot fully explain the observed vertical velocities, especially in the South-East Greenland, which confirms the need to refine our understanding of GIA contributions in this region.

SummaryIn this study, we adopt a horizontally layered model consisting of air, seawater and undersea porous rock and develop an analytically-based method to calculate the seismic and EM fields generated by undersea earthquakes. We conduct numerical simulations to investigate the characteristics of the EM response at the receivers located at the seafloor, in the seawater near the sea surface and in the air, respectively. The results show that two kinds of EM signals can be identified in the EM records at these receivers, namely, the early EM wave (seismic-to-EM conversion at the seafloor interface) arriving before the seismic waves and the coseismic EM fields with apparent speeds of the seismic waves. The EM signals observed at the seafloor are mostly stronger than those observed in the seawater and air near the sea surface. The method is applied to simulating the EM response to the 2022 Mw 7.3 earthquake that took place in the sea near Fukushima, Japan. At a receiver with 76 km epicentral distance at the seafloor, the predicted coseismic electric and magnetic signals reach 2 μV/m and 2 nT, respectively, which are within the detectability of the current EM equipment. This suggest a possibility to monitor the EM disturbances associated with undersea earthquakes and use them to serve the earthquake early warning, helping to mitigate the societal impact of large earthquakes.

SummaryLow-frequency laboratory measurements provide direct access to the elastic properties of samples within the seismic frequency band, offering calibration data for seismic survey analysis. Additionally, µCT imaging can quantify actual saturations and provides insights into phase distributions at the pore scale. To conduct laboratory triaxial measurements at seismic frequencies while simultaneously imaging the rock interior, we developed an X-ray transparent low-frequency apparatus. Our apparatus determines rock mechanical properties at seismic frequencies (0.5–150 Hz) and strain amplitudes (10−7–10−5), measuring Young’s modulus, Poisson’s ratio, and attenuation. In addition P- and S-wave velocities at ultrasonic frequencies are measured. We conducted imbibition-drainage experiments to assess the effect of saturation and patch size on seismic and ultrasonic elastic properties in sandstone. Additional tests with liquid and gaseous CO2 reveal the impact of partial CO2-gas saturation. The imbibition-drainage experiment demonstrated that P-wave velocity at ultrasonic frequencies was elevated during drainage and reduced during imbibition. Drainage caused patchy saturation, while imbibition resulted in uniform saturation. This implies that ultrasonic measurements, with wavelengths comparable to the pore fluid patch size, are likely influenced by scattering. In contrast, low-frequency measurements, where the wavelength surpasses the patch size, capture effective medium properties and therefore are not affected by scattering effects. The results of the CO2 test suggest that low-frequency measurements can detect even low gas saturations (4% gaseous CO2). In contrast, ultrasonic velocity measurements primarily reflect the response of the fully saturated sample at low gas saturations and do not indicate a reduction in velocity. Identifying fluid-solid interactions and estimating saturation via µCT imaging is crucial, especially with minimal gas presence. Our combined approach allows precise determination of elastic properties at seismic frequencies and shows the importance of low-frequency over ultrasonic measurements.

SummaryIntraplate earthquakes in stable continental regions exhibit diverse characteristics in terms of timing, spatial distribution, and magnitude. They are often unexpected, and their underlying physical mechanisms are not well understood. This complexity is particularly apparent in Norway, where seismicity is mostly localised on the continental margin and coastal areas. Various studies have attempted to explain the causes of seismicity in Norway by invoking different sources of stress, ranging from regional stress due to ridge push to local effects such as topography or deglaciation. In this study, we revisit these questions by investigating the distribution of seismicity in southwestern Norway using an enhanced earthquake catalogue. To achieve this, we revised the Norwegian National Seismic Network seismic catalogue from 2000 to 2023 and built a new catalogue using machine-learning-based techniques on data from a temporary seismic deployment in the region. Thanks to the increased station density during this deployment, we were also able to calculate new fault plane solutions that consistently showed a WNW-ESE direction for the most compressive axis. Furthermore, we demonstrate that seismicity in southwestern Norway, while diffuse, tends to be localised around the major crustal shear zones of the region, such as the Bergen Arc Shear Zone and the Hardangerfjord Shear Zone.

SummaryHaematite-bearing red beds are widespread across the Earth and play a pivotal role in palaeomagnetic studies. However, chemical remanent magnetisation (CRM) typically associated with authigenic haematite is not fully understood, which precludes more accurate interpretations of natural remanent magnetisation (NRM) in red beds. Here, we use electron microscopy, rock magnetism, and palaeomagnetism to investigate authigenic haematite in Early Triassic red beds in North China. Our findings reveal that the biotite-hosted haematite grains with grain sizes of several to tens of microns carry a significant portion of the NRM in these sedimentary rocks. We propose that these authigenic haematite particles primarily form during the early stages of diagenesis process. This authigenic haematite's growth is controlled by the crystal structure of the host biotite. Furthermore, this authigenic haematite displays high coercivity (> 100 mT) and high unblocking temperature (> 650 ° C), comparable to that of typical detrital haematite (30–1000 mT, > 650 ° C), which is usually the primary carrier of detrital remanent magnetisation (DRM) in such red beds. This study highlights the significance of combining mineralogical analysis with rock magnetism and palaeomagnetism to differentiate between CRM and DRM and thereby identify the primary NRM component within red beds. We hypothesize that the abundant iron supplied by biotite promotes the growth of authigenic haematite. This study illustrates the need to use caution when studying sedimentary NRM, particularly in rocks from source areas containing acidic igneous and metamorphic rocks (e.g. granite, diorite, and biotite gneiss) that contain a large proportion of iron-bearing minerals, such as the biotite observed in this study.