Geophysical Journal International

SummaryMagma transport in dikes is usually modelled by means of lubrication theory, assuming that magma properties are uniform across the dike. We explore the influence of cross-dike temperature heterogeneity on the dynamics of dike propagation using a quasi-2D model, derived from a full 2D model with an assumption of small width to length ratio. The model couples elastic fracture mechanics with multiphase magma flow, solving the governing equations using a hybrid numerical approach that combines the Displacement Discontinuity Method for elasticity with finite volume discretization for fluid flow and heat transfer. The model includes heat exchange with wall rocks, shear heating and latent heat release. It accounts for non-equilibrium magma crystallization, implementing temperature-dependent crystallization kinetics using an Arrhenius formulation for the relaxation timescale. As a case study, we simulate the ascent of a volatile-rich dacite from a source at 30 km depth. The distribution of temperature, crystallinity, and, thus, viscosity across the dike leads to a plug-like velocity profile with magma stagnation near the walls, substantially different from the parabolic Poiseuille flow assumed in classical lubrication theory. With temperature-dependent crystallization rate, rapid cooling of magma near the dike walls can generate a glassy chilled margin. The adjacent magma has higher crystallinity due to intermediate cooling rates, while the hotter core remains depleted in crystals throughout dike propagation. The dike propagates further and is thinner than predicted by (1D) lubrication theory because the low-viscosity core continues to facilitate vertical transport while the wall zones become progressively more viscous due to cooling and crystallization. The latent heat of crystallization can have a substantial impact in slowing down cooling and prolonging propagation. Other important factors include the characteristic crystal growth time, initial magma temperature and water content. Our quasi-2D approach bridges the gap between oversimplified 1D models and computationally expensive 3D simulations, providing a practical framework for investigating magma transport in silicic dikes.

SummaryThe 1957 Andreanof, Aleutian Is., earthquake (1957 March 9, 51.53°N, 175.63°W, d=25 km) is among the most enigmatic great earthquakes instrumentally recorded. The length of the aftershock area is very long (about 1,200 km), and tsunami excitation has been recently confirmed to be very extensive, yet its instrumental seismic magnitude Ms is only about 8.1 to 8.3. Detailed analyses of long-period surface waves in the past gave an Mw=8.4, and the seismic-tsunami disparity remains unresolved. The main difficulty in seismic studies is the absence of high-quality seismic data. Here we investigate the cause of this disparity by carefully analyzing some historical seismograms with modern digitization methods. We also take advantage of the 1996 Aleutian Is. earthquake (Mw=7.87) that occurred very close to the 1957 event. For the 1996 event, high-quality modern broad-band seismograms are available which can be effectively used as empirical Green’s functions for the analysis of the 1957 event. Using the Wiechert (Strasbourg, Uppsala), Milne-Shaw (Wellington), and Benioff (Uppsala, Pasadena) seismograms of the 1957 event, we could determine that the 1957 event had significant secondary excitation of long-period (150 s) waves during about 1,000 s following the first event. The Mw of the combined source is approximately 8.4. Because of the limited bandwidth of the old instruments, we cannot detect long-period energy beyond 150 s. However, the unusually long-lasting excitation over nearly 1,000 s suggests that the event had significant excitation at periods longer than 150 s with a much larger Mw for the total event. Although we cannot address this question quantitatively because of our band-limited data, our numerical experiment using a source with a slow component shows that if the time scale of the slow source is longer than 500 s, our data can be made compatible with an Mw =8.8 to 8.9 event, thereby reconciling the results from seismic and tsunami data.

AbstractThe rheological properties of the mantle govern plate tectonics and mantle convection, yet constraining the rheological parameters remains a significant challenge. Laboratory experiments are usually performed under different temperature-pressure-strain-rate conditions than those of natural environments, leading to substantial uncertainties when extrapolating the parameters to real-world conditions. While traditional Bayesian inversion with Monte Carlo sampling methods offers sufficient exploration of the parameter space and accurate inversion results, the excessive computational cost limits its practical application to complex nonlinear problems. To address these limitations, we integrate finite-difference-based geodynamic forward modeling with Automatic Differentiation (AD) to build a framework to invert non-linear rheological parameters. By incorporating multisource observational data, including surface velocities and topography, we are able to invert critical rheological parameters of the lithosphere and mantle, including the viscosity pre-exponential factor, activation energy, stress exponent, yield stress, and plate-interface viscosity. To validate the method, a series of models with different levels of complexity from single- to multiple-subduction systems and consideration of data noises are designed to generate synthetic data that are further used for inversion. Our method can successfully restore the rheological parameters under various conditions, with minimal errors between predicted and true values, underscoring its stability and broad applicability. In general, this study introduces a highly efficient and practical geodynamic forward and inverse modeling approach that can be used to infer the rheology of the mantle.

AbstractSeismic waveforms are essential for deciphering the subsurface structures of the Earth. Traditional methods for seismic waveform selection rely heavily on manual identification by experienced seismologists, which can be inconsistent and challenging when complex structures or huge amount of seismic data volumes are involved. Recent advancements in machine learning, particularly supervised learning techniques, have shown promising progress in addressing these challenges; however, their dependence on large labeled datasets limits their application to weak or rare seismic phases. In this study, we propose a new strategy using hierarchical clustering for seismic waveform identification, which does not need labeled dataset and minimizes extensive parameter settings. Our strategy is especially powerful when dealing with multiple waveform phases that may shift according to epicentral distance or may be distorted due to attenuation or other factors. We apply our strategy to identify various seismic wave of both P and S phases, especially those sampling deep Earth such as SKS-SKPdS and ScP phases. The results show that the strategy performs excellently and can identify different anomalous signals. Our approach empowers researchers to conduct more detailed studies in previously overlooked regions or datasets, thereby leading to a better understanding of deep Earth’s structures.

SUMMARYRock glaciers are specific landforms consisting of a mixture of rock debris, ice, liquid water, and air. In the Alps, active rock glaciers are generally found at high elevations above 2500 m. Active rock glaciers creep and can develop anomalous slide-like behaviors called destabilization. Induced polarization is a non-intrusive geophysical method that has proven to be sensitive to the hydrogeological properties of porous media. In August 2023, we performed four induced polarization profiles at Plan-du-Lac (Vanoise, France), on a multi-unit rock glacier complex with a front located at a low altitude of 2200 m). Our goal was to determine its architecture and its water and ice contents in relation with its activity rate. The survey included two transverse high-resolution profiles with a 5 m spacing between the electrodes and two other longitudinal profiles with a 20 m spacing between the electrodes allowing a depth of investigation of roughly 200 m to image the rock glacier from its terminal front up to its root. The conductivity and normalized chargeability tomograms were inverted and then used to get the water content and cation exchange capacity (a proxy for the clay content) tomograms. In most of the units, ice has disappeared and the landforms associated with the former rock glacier were characterized by low water and clay contents with respect to the basement. This was consistent with these units being mostly formed by rock debris with a low water saturation except at their bases, which are water-saturated. Ice remains were found at the roots of the rock glacier, with a volume content up to ∼10 per cent (vol. per cent) for profile P2 and 16 per cent for profile P4. The roots of the rock glacier complex were still creeping as shown by InSAR data. This case study demonstrates the usefulness of the induced polarization method to quantitatively characterize gravitational instabilities associated with coarse materials and transitional rock glaciers.

SUMMARYIn October 2023, an earthquake sequence comprising four ∼ MW 6.0 events struck Herat Province in northwestern Afghanistan, causing severe casualties and property losses. The geometry of seismogenic faults and the mechanisms of the earthquake sequence are essential for regional seismic hazard assessment, but still remain poorly constrained. With Interferometric Synthetic Aperture Radar (InSAR) techniques, we extracted high-resolution co- and early post-seismic deformations of the events. Through a two-step inversion method, we inferred the geometry of the causative faults and the distributed slip models. The earthquake sequence ruptured two intersecting low-dip thrust faults, indicating that the complex geometry may have played a key role in controlling the propagation of the events. The ruptures of the four major events are clearly imaged at depths of 1-10 km without reaching the surface, showing a pattern of first spreading westwards, then jumping eastwards to the bend segment, and finally rupturing an adjacent fault. Post-seismic deformation further reveals reactivation of a secondary fault splay which underwent afterslip. Shallow afterslip up-dip of the co-seismic rupture dominates post-seismic deformation during 10 months following the earthquake sequence. Relying on the evolution of afterslip, we infer that significant rate-strengthening property in the shallow bend section may have hindered further co-seismic rupture propagation. Combining obtained results and the complex geological setting of the Herat region, we suggest that the earthquake sequence reflects N-S crustal shortening between two branches of the western Herat Fault System.

SUMMARYIn the last decade, the Dynamic Stern Layer (DSL) model has proven to be a reliable petrophysical model to comprehend induced polarization data at various scales from the representative elementary volume of a porous rock to the interpretation of field data. Preliminary works have demonstrated that such model can be extended to understand the induced polarization properties of ice-bearing rocks and to interpret field-acquired induced polarization data in the context of permafrost. That being said, the direct effect of ice was let aside. We first review the DSL model in presence of ice and discuss the role of ice as an interfacial protonic dirty semi-conductor in the complex conductivity spectra with an emphasis on the role of the complex-valued surface conductivity of ice crystals above 1 Hz. We propose a new combined polarization model including indirect and direct ice effects. By direct effects, we mean the effects associated with change in the liquid water content and salinity of the pore water. By direct effect, we mean the role of the interfacial properties of the ice surface and liquid water is still present in the pore space of the porous composite. In this case, the electrical current is not expected to cross the ice crystals. Instead, it would polarize the surface of the ice crystals and generate a very high chargeability that can reach one depending on the value of the volumetric content of ice. We apply the DSL model to a new set of complex conductivity spectra obtained in the frequency range 10 mHz-45 kHz using a collection of 25 rock samples including metamorphic and sedimentary rocks in the temperature range + 15/+20°C to -10/-15°C. We observe that the model explains very well the observed data in the low-frequency range (10 mHz-1 Hz) without any direct contribution of ice. In the high frequency range (above 1 Hz), we observe a weak contribution possibly associated with the contribution of ice crystals. We establish under what conditions the direct contribution of ice can be neglected. We also investigate the role of porosity, cation exchange capacity, and freezing curve parameters on the complex conductivity spectra of crystalline and non-crystalline rocks during freezing. Laboratory experiments demonstrate that in most field conditions including permafrost conditions, surface conductivity associated with conduction on the surface of clay minerals (and alumino-silicates in general) is expected to dominate the overall conductivity response. Therefore Archie’s law cannot be used as a conductivity equation in this context because of the contribution of surface conductivity. A large experimental and field dataset at the Aiguille du Midi (3842 m a.s.l., French Alps) for the resistivity versus temperature data of granitic rocks demonstrates the role of surface conductivity in the overall conductivity of the rock.

SummaryGangwon Province, located in the central part of the Korean Peninsula, features northeast–southwest faults and tectonic structures formed by plutonic intrusions. Despite decades of geological investigations from near-surface to the upper crust in Gangwon Province, the lithospheric structure of this region remains poorly understood. The primary objective of this study is to identify velocity anomalies potentially associated with plutonic intrusions and to elucidate the formation processes and mechanisms governing the crustal and upper mantle structures in this region. We employed Helmholtz tomography to generate phase-velocity maps for periods of 10–40 s using a dense seismic network of 101 stations. These maps were subsequently inverted to obtain an S-wave velocity model from the upper crust to the uppermost mantle. Our results reveal northeast–southwest-trending low-velocity anomalies along major faults in central to northern Gangwon Province (i.e. eastern Gyeonggi Massif), extending to depths of approximately 25–30 km. These low-velocity anomalies align with the orientations of Jurassic granitoid intrusions formed through partial melting processes. Additionally, we identified other low-velocity anomalies, likely formed by Late Cretaceous intrusions, which are oriented perpendicular to the major faults. In contrast, the southeastern Gangwon Province (i.e. Taebaeksan Basin) exhibits a distinctly different velocity structure, lacking features indicative of granitic intrusions and showing low-velocity anomalies confined to shallow depths. The pronounced low-velocity anomalies observed at depths of 5–10 km in Taebaeksan Basin are attributed to a complex fault zone influenced by Permo-Triassic collisional orogeny.

SummaryWith the deployment of dense linear seismic arrays, teleseismic waves are playing an increasingly important role in studying deep structures beneath seismic stations. However, despite significant advancements in high-performance computing, simulating high-frequency teleseismic waves (above 1 Hz) in a 3D model on a global scale remains challenging. To address this issue, previous studies have developed hybrid methods that utilize the displacement representation theorem to equivalently transform stresses and velocities simulated in a 1D or 3D global reference earth model into body force and moment rate density tensor sources for input into a 3D region model. Although previous hybrid methods have incorporated the free surface, the treatment of two types of equivalent sources at this boundary, particularly the equivalent moment rate density sources, has not been fully addressed. Neglecting the influence of the free surface condition and directly adding the equivalent sources at the free surface may lead to spurious waves. To resolve this, we develop a hybrid simulation method considering the effects of the free surface condition. By setting the relevant components of the equivalent sources on the free surface to zero, the method effectively reduces artifacts caused by coupling effects. We then propose the QSSP-CGFD3D hybrid method, which includes this free surface boundary correction, for simulating teleseismic waves in a 3D receiver-side model. We validate the accuracy and effectiveness of the hybrid method for calculating P-waves, S-waves, and surface waves in the AK135 model. We also apply the method to the fault zone region, where the results show that the fault zone causes arrival time delays and amplitude amplifications of teleseismic P-waves. These effects can be used to infer structural parameters of fault zones. Furthermore, we employ the QSSP-CGFD3D hybrid method to simulate the influence of undulated interface within the crustal structure on telesesimic waveforms, demonstrating its potential for receiver function analysis. The proposed hybrid method demonstrates significant potential for studying structures beneath seismic arrays, and holds promise for advancing our understanding of such features.

SummaryBack-projection inversions of teleseismic waveforms for the images of rupture progression in great earthquakes have become a popular tool in earthquake studies. However, verifying the trustworthiness of the obtained images in synthetic tests, in which forward-problem data from known propagating ruptures on finite faults are inverted and compared with the true images, have been disproportionally lacking. Such validations for known rupture geometries in a homogeneous medium provide the best-scenario probe into the theoretical ability of the method to resolve the true faulting kinematics. Unambiguous identification of the true source of radiation is possible if alternative trial subfaults, not representing the real emitter, shift the wave-arrival time by different amounts at different stations, resulting in no significant stack achieved for them. This condition is quantitatively expressed in the value of the dimensionless uniqueness coefficient q, which must exceed unity for the inversion to become unique. The criterion is not satisfied for the faults of finite dimensions, precluding reliable determination of rupture progression and speed for them, no matter how wide the coverage of the azimuths from the fault to the stations and how many stations in the network. Unambiguous determination of the exact location is possible for point-source radiators with widening of the azimuthal coverage: the reduction in ambiguity is seen in the progressive improvement in the correctness of the images as the uniqueness coefficient increases to the values greater than one. Fictitious moving, linearly aligned sources appear, increasing in number, as the value of the coefficient gradually drops.

SummaryOlivine, orthopyroxene and clinopyroxene are the most common minerals in ultramafic rocks, their modal contents and crystallographic preferred orientations (CPOs) are main factors determine the total seismic properties (e.g. seismic velocity and anisotropy). Red Hills Massif is the main part of Dun Mountain Ophiolite Belt, and includes all six most typical olivine CPO types (A-E & AG-type) and various ultramafic rock types. In this case study, 6 paired harzburgite and dunite (Series 1 samples) and olivine/clinopyroxene-related ultramafic rocks (3 dunites, 2 wehrlites, 2 olivine clinopyroxenites and 1 clinopyroxenite) are selected to proceed detailed seismic characteristics analysis. Seismic distributions of Series 1 and 2 peridotites are based on their olivine CPO characteristics while the distributions of Series 2 (olivine) clinopyroxenites are based on their clinopyroxene CPO characteristics. Due to the addition of orthopyroxene, almost all harzburgites have slower P-wave velocity and less P/S-wave anisotropy than dunites in the same pair. On the other hand, both olivine & clinopyroxene CPO combinations and modal contents of Series 2 ultramafic rocks are various. With increasing modal proportions of orthopyroxene or clinopyroxene, seismic properties are not necessarily decreasing. When the [001]OL crystallographic axis in harzburgite develops a girdle fabric, or when the [100]OL and [001]CPX orientations exhibit misalignment with respect to the lineation direction, a significant reduction in seismic wave velocities is observed. Concurrently, seismic anisotropy magnitudes demonstrate marked enhancement under orthopyroxene- or clinopyroxene-dominant conditions (>60 per cent modal abundance). The intricate olivine and ortho-/clinopyroxene CPO patterns and related seismic characteristics preserved beneath the Red Hills Massif are thought to pre-date initiation of the Alpine Fault (∼25 Ma). This interpretation is supported by the similarity in olivine grain size between the Red Hills peridotites and the lithospheric mantle beneath West Otago, implying a shared pre-25 Ma mantle domain. Prior to Alpine Fault offset, these two regions were adjacent, and their lithospheric high-velocity seismic signatures remain correspondingly alike. Once the initiation of Alpine Fault kinematics started, the observed anomalous SKS azimuth proximal to the study area may reflect either: (1) the development of protomylonitic textures analogous to those documented in the West Otago lithospheric mantle, or (2) the preservation of pre-existing deformation fabrics. Conversely, lithospheric fast shear wave splitting directions, when combined with gravity data and rock density constraints, can potentially resolve the dominant olivine CPO type(s), the relative proportions of olivine and ortho-/clinopyroxene, and their combined CPO patterns.

AbstractTomography based on full waveforms is an important tool for characterising the subsurface. However, systemic artefacts in the recorded data must be removed prior to the imaging process to fully utilise the information contained in the data. Especially for near-surface surveys, the coupling between sources/receivers and the medium can introduce significant distortions in the recorded data. We present two novel time-domain FWI approaches that account for the interaction between sources, receivers, and the subsurface. The first approach does not impose any restrictions on the shape of the source wavelet, except that it must be compact in time (TC inversion) such that the computation of synthetic seismograms is feasible. In the second approach, we assume that the source wavelet can be approximated with a Ricker wavelet, and we invert only for the three parameters describing a Ricker wavelet (SP inversion). Both algorithms have been tested with synthetic crosshole data. The SP approach is slightly superior to the TC inversion when the true wavelet is indeed a Ricker wavelet. However, the TC inversion outperforms the SP approach when the true source wavelet is not well approximated with a Ricker wavelet. This is demonstrated with a field data set acquired in boreholes at a CO2injection test site in Svelvik (Norway).

SummaryThe probability that any given time interval of duration τ is occupied by one earthquake, or more, characterizes the temporal statistics of seismic sequences and, therefore, the temporal clustering of events. The occupation probability reveals the fractal behaviour of seismic sequences, $\Phi (\tau ) \sim \tau ^{1 - D_\tau }$, at short times, defining a temporal fractal dimension, Dτ, for seismic events. We introduce an empirical model of the occupation probability, parametrized by the fractal dimension and two other parameters. We use the mathematical relationship between the occupation probability and the inter-event time probability density to develop intuition about the model parameters and to compare the proposed model with the conventional gamma model for inter-event times. Our model captures the statistical properties of a wide range of seismic sequences where the gamma model fails, from slow-slip driven swarms to burst-like episodes at the bottom of the seismogenic zone and low-frequency earthquakes. Using real and synthetic catalogues, we find that the model parameters are related to the degree of intermittency of the seismic activity, the characteristics of the bursts and the time scale of the quiescent periods. Measuring and modelling the occupation probability constitutes a valuable tool to categorize seismic sequences over a wide spectrum of seismic occurrence patterns. When applied to the new generation of earthquake catalogues, this empirical method highlights intermittency and fractal bursts, challenging conventional seismicity models and emphasizing the need to refine them to better capture these key features.

SummaryNowadays, many joint inversions are carried out to understand the earthquake source process. In the joint inversion analysis, we have to determine the relative weights among different datasets in addition to the regularization term, such as smoothing. Akaike’s Bayesian Information Criterion (ABIC) is known to be useful to find the appropriate values of such hyperparameters. This study proposes a method to jointly invert tsunami, GNSS, and SAR data using ABIC to construct a finite fault model. We demonstrate our inversion scheme in the case of the 2024 Noto Peninsula earthquake, whose fault geometry is still under discussion. Since the dip angle of the fault can also be considered as a hyperparameter, we evaluate three types of dip angles and estimate an appropriate value based on ABIC. In other words, our inversion scheme utilizes ABIC to determine the dip angle, the weights among datasets, and the spatial smoothness of fault slip. Our fault model indicates that (1) listric fault, varying the dip angle with depth, is the most appropriate among the ones we proposed, (2) the largest slip is on the fault under the northwestern corner of the peninsula, and (3) coseismic fault slip extends to offshore faults east of the peninsula. In the case of the listric fault, ABIC values GNSS and SAR data, which improves the agreement of the on-land coseismic displacement while also reproducing tsunami data. We also find that analyzing tsunami records in the frequency domain helps to obtain a robust inversion result when employing ABIC.

SummaryThe elastic deformation of Earth’s surface caused by internal mass distribution varies significantly across loading models, especially in high-precision applications. Although several studies have applied loading corrections to Global Navigation Satellite System (GNSS) time series in mainland China, discrepancies between models, particularly those involving Gravity Recovery and Climate Experiment (GRACE) data and its downscaled derivatives, remain insufficiently explored. Moreover, previous research has not comprehensively assessed vertical crustal deformation after applying different environmental loading corrections. This study systematically evaluates the correction effects of various environmental loading models on GNSS vertical displacements across mainland China, generating vertical velocity maps along with their associated uncertainties. The results show that hydrological loading (HYDL) has the most significant impact on GNSS vertical displacements, whereas non-tidal oceanic loading (NTOL) has the least effect. Substantial differences exist between various HYDL models, while discrepancies between non-tidal atmospheric loading (NTAL) and NTOL models are relatively minor. A comparison of correction effects between the HYDL model and GRACE data reveals that the HYDL model offers more accurate corrections, whereas downscaled GRACE data demonstrates improved performance, underscoring its potential advantages. After applying loading corrections and filtering common mode error (CME), the uncertainty in GNSS vertical velocity is notably reduced, although velocity variation remains small. This effect is also evident in seasonal variations. Furthermore, a comparison of vertical land motion (VLM) constrained by different HYDL models and downscaled GRACE data with VLM constrained by an independent land-water fusion model reveals higher consistency between the downscaled GRACE data and the independent model in the North China and northwestern Tianshan regions, suggesting that VLM derived from downscaled GRACE data may be more reliable. We also quantify the combined impact of geocenter motion and glacial isostatic adjustment (GIA) on the VLM trend across mainland China, estimated it at approximately 0.13 mm/yr. While the spatial characteristics of the VLM trend show minimal changes after correction, its intensity is significantly affected. This study provides crucial insights into the correcting of environmental loading effects in GNSS vertical displacements and contributes the latest observational results on vertical crustal deformation in mainland China.

SummaryThis study presents an enhanced method for integrating Global Navigation Satellite System (GNSS) and Interferometric Synthetic Aperture Radar (InSAR), referred to as Iterative Decomposition-based GNSS-enhanced InSAR (ID-GInSAR), to address both spatially correlated components (SCCs) and topographically correlated components (TCCs) in interferogram errors. While traditional GInSAR (GNSS-enhanced InSAR) is effective in mitigating long-wavelength SCCs, it often overlooks TCCs, which are particularly significant in regions with steep topographic gradients. The proposed ID-GInSAR approach employs an iterative decomposition process to decouple and independently model SCCs and TCCs, utilizing a combination of exponential and statistical models. The method is validated using Sentinel-1 SAR and GNSS data from California’s southern Central Valley. Results demonstrate that ID-GInSAR significantly lowers noise in interferograms, enhances the robustness of displacement time series, and improves the accuracy of co-seismic deformation and velocity field estimates. Specifically, ID-GInSAR reduces the root mean square (RMS) between GNSS and InSAR by up to 55% in individual interferograms and by an average of 30.4% in displacement time series compared to traditional GInSAR methods. Furthermore, ID-GInSAR effectively highlights subtle transient deformation, such as coseismic offsets, and provides more robust velocity fields over shorter time spans (less than three years). Finally, we compare our method with other approaches, including Remove/Filter/Restore (RFR) and GACOS, and discuss their applicability scenarios. Collectively, ID-GInSAR provides an alternative integration method for regions with complex topography where ground-based GNSS observations are available.

SummaryRupture directivity significantly increases horizontal peak ground acceleration, elongates aftershock clouds, and enlarges meizoseismal areas beyond the fault end in front of the direction of rupture propagation. In this study, we examine the directivity of 25 moderate to large earthquakes (Mw ≥ 6) from 1968 to 2017 in the Iranian plateau by employing relocated earthquake clusters, mapped surface ruptures, focal mechanisms of earthquakes, slip distribution models, spatial distribution of Peak Ground Acceleration (PGA) amplitudes and macroseismic effects. The methodology overcomes the lack of dense seismic networks required to study directivity using methods based on the azimuthal variation of the spectrum of seismic waves. We show that 16 out of the 25 (i.e., 64%) of the earthquakes investigated have mostly unidirectional rupture. This implies that unidirectional ruptures in a slow deforming continental collision zone such as the Iranian Plateau is only slightly less common than those observed globally. With the understanding that unidirectional rupture increases the probability of ground shaking off the termination of the causative faults, our findings highlight the importance of considering the directivity effect in earthquake hazard assessment in Iran and also in other slow deforming continental regions.

SummaryEvidence from seismic studies, mineral physics, thermal evolution models and geomagnetic observations is inconclusive about the presence of a stably stratified layer at the top of the Earth’s fluid outer core. Such a convectively stable layer could have a strong influence on the internal fluid waves propagating underneath the core-mantle boundary (CMB) that are used to probe the outermost region of the core through the wave interaction with the geomagnetic field and the rotation of the mantle. Here, we numerically investigate the effect of a top stable layer on the outer core fluid waves by calculating the eigenmodes in a neutrally stratified sphere permeated by a magnetic field with and without a top stable layer. We use a numerical model, assuming a flow with an m-fold azimuthal symmetry, that allows for radial motions across the lower boundary of the stable layer and angular momentum exchanges across the CMB through viscous and electromagnetic coupling. On interannual time-scales, we find torsional Alfvén waves that are only marginally affected by weak to moderate stratification strength in the outer layer. At decadal time-scales similarly weak stable layers promote the appearance of waves that propagate primarily within the stable layer itself and resemble Magneto-Archimedes-Coriolis (MAC) waves, even though they interact with the adiabatic fluid core below. These waves can exert viscous and electromagnetic torques on the mantle that are several orders of magnitude larger than those in the neutrally stratified case.

SummaryA novel time-lapse modelling scheme for Airborne Electromagnetics (AEM) monitoring datasets is presented, using data from multiple surveys applied to study the hydro-related evolution of the Bookpurnong floodplain in South Australia. Additionally, it introduces a new wide-ranging approach for this type of study, incorporating new processing, validation, and interpretation tools.Time-Lapse studies are widespread in the literature but are not commonly applied to model EM data, particularly AEM data. This is linked to the challenges of performing overlapping data acquisition with inductive systems. The key features of the new time-lapse scheme, which address these issues, include the definition of independent forward and model meshes, essential for considering discrepancies in the location of soundings which arise in multitemporal AEM data acquisition, and the incorporation of system flight height in the inversion. This proved crucial for achieving satisfactory data fitting and limiting artifact propagation in the time-lapse models.Additionally, a novel processing workflow for AEM multitemporal datasets is presented. This has proven important for effectively processing the multitemporal datasets, which presents new challenges in identifying noise coupling arising from the use of different systems across vintages of data, possible variations in acquisition settings operated by different field crews, and changes in subsurface resistivity in the survey area. Results generated from the time-lapse modelling are evaluated with an Independent Hydrogeological Validation (IHV), designed to support the geophysical models validation and interpretation by providing a first-step hydrogeological evaluation.At Bookpurnong, along a sector of the Murray River floodplain, multitemporal AEM surveys were collected in 2015, 2022 and 2024, to study natural and engineered changes in the groundwater system over time. The time-lapse models show significantly smaller variations compared to those determined with individually modelled survey data sets, while delineating sharply bounded changes in resistivity across the floodplain. This demonstrates the effectiveness of the new time-lapse scheme in minimizing inversion variations typically encountered with independently modelled results affected by larger equivalence issues.Here, AEM models are first compared with resistivity borehole measurements, revealing a strong match between the two methodologies and spatial variations in resistivity consistent with a meandering river across the floodplain. These variations are further validated and interpreted using the IHV approach, which revealed a direct correlation between the hydrological stress of the Murray River and the response of shallow aquifers. Additionally, time-lapse geophysical models, combined with a hydrostratigraphic analysis, allow for a direct correlation between shallow and deep hydrogeological responses.We believe that the time-lapse methodology described here can be widely applied to multitemporal studies using AEM datasets, enabling the study of a broad range of natural processes with great accuracy and at the basin scale.

SummaryThe complex three-dimensional (3D) geometry of active faults plays a crucial role in controlling earthquake location, extent, and rupture behavior, making the accurate representation of fault models essential. Fault structures are typically interpreted manually from relocated hypocenters and interpolated to generate 3D fault surfaces. However, this process is often non-unique and uncertain due to the uneven spatial distribution of earthquake hypocenters, the subjectivity of manual interpretation, and the complexity of non-planar faults. To address these challenges, we developed a method that combines adaptive threshold hierarchical clustering with quantitative evaluation to automatically and effectively construct 3D models of seismogenic faults. This method utilizes the nearest neighbor index (NNI) to determine whether seismic activity exhibits clustering characteristics indicative of fault structures. Adaptive threshold hierarchical clustering is subsequently applied to identify small earthquake clusters associated with each fault. High-density 3D automatic slicing ensures robust fitting of fault lines, and in combination with surface rupture data, discrete smooth interpolation (DSI) is used to construct a 3D fault model. For each fault, we calculate distances from small earthquake clusters to the 3D fault structure and analyze their spatial distribution using kernel density estimation (KDE) to optimize the model for a near-symmetric distribution of small earthquake clusters on both sides of the fault. We applied this method to the 2013 Ms 7.0 and 2022 Ms 6.1 earthquakes in southern Longmenshan, Sichuan, China, refining the 3D seismogenic fault models for both events. Additionally, we constructed a 3D fault model for the 2019 Ridgecrest Mw 7.1 earthquake sequence using the same approach. The results indicate that this method is applicable to both individual faults and multiple intersecting fault systems. Compared to traditional manual modeling approaches, our method significantly enhances the identification of small earthquake clusters, reduces reliance on manual interpretation, increases modeling efficiency, and minimizes errors. This innovative modeling technique advances the 3D geometric construction of complex active faults and is adaptable to a wide range of seismic research applications.