Geophysical Journal International

SummaryGlobal Navigation Satellite System (GNSS) plays a fundamental role in monitoring time-dependent ground displacement. However, GNSS daily position time series can often contain significant outliers, reaching up to several centimeters. These are likely of non-tectonic origin, and, if not properly accounted for, they can significantly impact the accuracy and dependability of the estimation of key parameters for geophysical analyses, such as long-term velocities and transient deformations. Characterizing these outliers can provide information about their possible sources and help us implementing mitigation strategies. Asymmetric outliers, i.e., those characterized by a primary direction, therefore occurring on one side of the mean time series, are of particular interest since they could point to the presence of recurring or repeatable sources of error. Their key features and potential causes are, however, still not fully analyzed and understood. We analyze asymmetric outliers in thousands of GNSS time series across three regions – Central-Southern Italy, New Zealand, and the Western U.S. – using data from multiple processing centers, and we reveal some persistent features among all datasets. Tens of the analyzed sites exhibit hundreds of large outliers (10-50 mm), far exceeding typical position uncertainties (∼1-6 mm). Remarkably, the outliers are numerous in the horizontal component, and tend to occur near mountainous regions, with preferred direction roughly orthogonal to the local topography. The results consistency across different datasets and instrumental features suggest a physical origin for these outliers rather than a specific processing approach or instrumental configuration. Further analyses at local scales align with previous studies linking skewed position errors to uncorrected tropospheric delays driven by the coupling between atmospheric conditions and local terrain (e.g., trapped lee waves). However, other factors – such as multipath, snow accumulation on GNSS antennas or obstructed sky visibility – could also contribute to the observed asymmetric outliers. We explore mitigation strategies at both processing and post-processing stages, but further analyses and more sophisticated approaches, such as high-resolution tropospheric modeling, are needed to better understand the involved processes and achieve meaningful improvements.

SummaryThe Gravity Recovery and Climate Experiment (GRACE) and its Follow-On mission provide essential observations of Earth’s surface mass redistribution. However, inherent north-south striping noise in the GRACE spherical harmonic (SH) products limits their application at sub-basin scales. To address this, we introduce a novel spatial domain decorrelation filter, the Physical-Informed Spatial Pattern (PISP) filter, which leverages the structured physical characteristics of the noise for its precise identification and removal. Comprehensive numerical experiments validated that PISP effectively eliminates striping noise globally and yields a consistent noise background across latitudes, with noise reduced to a uniform level in more than 90% of the months examined and with stable performance under strong-noise conditions. In a case study of water storage variations in Lake Victoria, PISP preserves the primary signal amplitude and reduces the root-mean-square error relative to reference data to 5.84 cm after spatial smoothing, outperforming the 6.81 cm achieved by the MVMDS + DDK6. Furthermore, for three earthquakes with magnitudes exceeding 8.8, PISP effectively removes striping noise using regional masking, successfully recovering the co-seismic signal morphology. By further verifying the method’s stability across various noise scenarios, the results demonstrate PISP’s potential for future global research integrating multi-satellite gravity data.

SummarySlowslip events (SSEs) modulate the earthquake cycle in subduction zones, yet understanding their physics remains challenging due to sparse observations and high computational cost of physicsbased simulations. We present a scientific machine-learning approach using a data-driven reduced-order modelling (ROM) framework to efficiently simulate the SSE cycle governed by rate-and-state friction in a Cascadia-like 2D subduction setting. Our approach projects fault slip, sliprate, and state-variable trajectories onto a splinebased latent space, which is subsequently emulated using properorthogonal decomposition and radialbasisfunction interpolation. Achieving a speedup of ∼360, 000 × compared to volumetric simulations, the ROMs enable comprehensive parameter exploration and Bayesian Markov chain Monte Carlo (MCMC) inversion. By smoothly interpolating between the physics-based simulations, the ROMs reveal complex dependencies that might be overlooked with coarser parameter space sampling. Our analysis reveals complex, non-linear dependencies of SSE characteristics on the width and magnitude of the deep, low-effective-normal-stress region while holding frictional parameters constant. We show that, to first order, the recurrence time of SSEs is linearly dependent on the normalized fault width, defined as the SSE zone width divided by the characteristic nucleation length, in agreement with previous work. However, at second order, the recurrence interval can change more rapidly with small variations in the SSE zone width. We identify a region of steep, non-linear dependence of the recurrence interval on the normalized fault width, which we attribute to the extent of the velocity-weakening portion of the subduction interface that produces SSEs. Our MCMC inversion constrained by Northern Cascadia SSEs observations indicates near-lithostatic pore fluid pressure (99.6 ± 0.17% lithostatic) and positions the upper frictional transition zone at 30.4 ± 2.8 km depth, consistent with geophysical observations. The inversion resolves the deep SSE-portion of the slab spanning 45 ± 16 km with low effective normal stress of 3.8 ± 1.4 MPa. We discuss how varying the fault frictional parameters affects the MCMC-inverted parameter values. This framework provides a new tool for advancing the physics-based understanding of SSEs and subduction zone faulting mechanics. By systematically linking megathrust properties such as fluid pressure and fault strength to rate-and-state friction governed slow slip cycle characteristics, such as recurrence interval, our approach helps to constrain the first and second-order physics-based controls and the uncertainties of how subduction zones slip.

SummaryDuring the last glacial period, the Black Sea developed oxic bottom waters that generally preserved detrital titano-/magnetite, yet greigite (Fe₃S₄) occurs sporadically within nodules in these otherwise oxic sediments. To clarify their origin and significance, we combined rock magnetism, X-ray diffraction, and scanning electron microscopy with energy-dispersive X-ray spectroscopy on nodules from core MSM33-55-1 and screened fifteen nearby cores using the ratio of saturation isothermal remanent magnetization to low-field magnetic susceptibility (SIRM/κLF). SIRM/κLF thresholds were applied to classify samples as greigite-bearing and greigite-free, then the greigite-bearing sample proportion was calculated in 0.3-kyr bins from 70 to 20 ka. The nodules have S-rich interiors (elemental sulphur plus greigite) coated by Fe-hydroxide rims (goethite), indicating formation in localized sulphidic microenvironments embedded in otherwise oxic sediment. Palaeomagnetic comparisons indicate that greigite-bearing samples track greigite-free (detrital titano-/magnetite) records for relative palaeointensity (RPI) and inclination (zero lag for inclination; ∼0.5 kyr lag for RPI), which implies that greigite formed shortly after deposition and acquired an early-diagenetic chemical remanent magnetization near-synchronous with a detrital post-depositional remanent magnetization. Greigite-bearing sample proportion increases during interstadials with modest but significant positive correlations with total organic carbon (TOC) content (and TOC/κLF) and a negative correlation with κLF, consistent with enhanced organic-matter supply relative to detrital input promoting sulphidic microenvironments. Together, these results demonstrate that greigite-bearing nodules can be incorporated into palaeomagnetic reconstructions when carefully screened, and that they serve as markers of micro-scale sulphidisation coupled to conditions that favour organic-carbon retention in glacial Black Sea sediments.

SummaryIn arid and semi-arid climates, it is critical to assess the state of the soil in terms of water content and water salinity. The use of geophysics, and geoelectrical methods in particular, to this end faces the challenge of discriminating between the effects of water content and pore water salinity on soil electrical resistivity, since these two factors are both inversely related to resistivity. The heterogeneity of soils, with its possible varying clay content, makes the interpretation even more complex. We have investigated the combined effect of water saturation, pore water salinity and clay content using Spectral Induced Polarization (SIP) on controlled laboratory samples. The experiments contribute to the development and application of the SIP method in the area of small-scale data acquisition and processing for hydrogeophysical and environmental purposes. In our experimental setup the above-mentioned three variables were gradually modified under controlled conditions. Ca-montmorillonite and very fine to coarse sand were mixed during multiple dry-wet mixing cycles in order to create artificial soil samples that mimic natural soils. In total, 4 samples were used with clay content varying from 2 to 8 mass per cent of clay. The other two variables were changed as well: water saturation ranging from 100 per cent to 10 per cent in 9–15 steps, and electrical conductivity of the pore water ranging from 0.05 to 0.7 S m−1 in 4 steps. The statistical analysis of the results indicates that there is a significant positive correlation between quadrature conductivity and the three variables. The obtained data was fitted using a double Cole-Cole model. The analysis of the obtained Cole-Cole parameters along the three variables shows promising results for the separation between the effects of pore water salinity and water content, thus paving the way for fruitful field applications in arid and semi-arid environments.

SummaryIn the context of detecting shallow, localized seismic velocity variations, we assess the sensitivity of wavefield gradients, specifically normal strain and rotation, relative to traditional seismological observables such as displacement, velocity, and acceleration. We begin with a simple single-scattering analysis of the displacement wavefield and its gradient, and introduce the Proximity Field Test (PFT) as a straightforward and effective tool for subsurface detection. We then perform 3D elastic simulations with SEM46, a spectral-element code modified to directly output strain and rotation, and analyze two case studies involving localized shallow velocity changes in a homogeneous medium. In the first case, P- and S-wave velocities and density are varied by 10 per cent relative to the background, whereas in the second case the local velocity decrease is 70 per cent. By comparing waveforms in the reference medium with those obtained after introducing the heterogeneity, we show that the joint analysis of displacement and gradient measurements enables efficient detection of subsurface anomalies. Finally, we validate the approach with a field experiment combining geophones and DAS measurements, aimed at detecting a buried concrete foundation associated with an approximate 70 per cent positive velocity contrast in the shallow subsurface. The anomaly is clearly identified with minimal processing, demonstrating the practical potential of the proposed method.

SummaryGeothermal energy production emits significant amounts of hydrogen sulfide (H2S). A strategy to mitigate the emissions is to reinject the H2S into basaltic formations, where it reacts with the rock to form pyrite. Due to the polarization properties of pyrite, the direct current resistivity (DC) and induced polarization (IP) geophysical method (i.e., DCIP) has shown the potential for monitoring H2S mineral storage. However, field applications of DCIP monitoring have been limited by the low spatial coverage of wireline logging and by the ambiguity in interpreting IP signals due to multiple processes that contribute to the polarization response. This study integrates DCIP with field-scale reactive transport modeling, utilizing both synthetic modeling and field investigations, to assess the ability of surface and cross-hole DCIP to monitor H2S mineral storage at the Nesjavellir study site in Iceland. Two surface DCIP datasets were collected at Nesjavellir, with six months of continuous H2S injection between them. Time-lapse inversions, performed using a novel gridding scheme that accounts for electrode misplacement between the two DCIP surveys, recover no significant IP changes beyond data noise. Interpreting these results alongside the reactive transport results finds that pyrite mineralization during the six-month injection period is too small and too deep to be resolved by surface DCIP time-lapse surveying, highlighting the benefit of a joint geophysical-geochemical interpretation approach. Conversely, joint DCIP-reactive transport synthetic modeling shows the potential of cross-hole DCIP for monitoring long-term H2S mineral storage, provided that data noise is low and sufficient H2S is injected. The reactive transport models also provide insight into the mechanisms contributing to the polarization response, demonstrating that pyrite mineralization is the primary contributor to the polarization response, with minimal contribution from other minerals such as smectites and iron oxides. However, existing petrophysical relationships are simplistic, which adds uncertainty to the interpretation of the DCIP signal and the quantification of pyrite mineralization. Additionally, smectite formation has been shown to decrease both the polarization signals and the quality of the IP data due to its electrically conductive properties. At Nesjavellir, a decrease in the DC resistivity from 1925 to 325 Ωm is observed, attributed to the disposal of warm wastewater. This decrease is identified by comparing resistivity data collected in this study to historical vertical electrical sounding data collected prior to geothermal development. Lastly, petrophysical relationships linking DC resistivity, smectite content, and permeability suggest a high basalt fracture permeability of 7.9× 10−12 m2, which agrees with values recovered through flow model calibrations. This result demonstrates the value of geophysical surveying not only for monitoring but also for constraining key parameters in reactive transport simulations.

SummaryMegathrust earthquakes occur at subduction zones and produce many of the largest magnitude earthquakes on record. The ruptures of megathrust earthquakes have classically been thought to nucleate at asperities, zones on the subduction interface of rheological strength which accumulate strain as the plates move. Imaging such asperities would critically inform seismic hazard analyses, to the benefit of millions who live adjacent to major subduction boundaries. We created a 3D shear wave attenuation (QS) tomography of the 2015 MW 8.3 Illapel, Chile, earthquake rupture region by measuring QS using a body wave S/P spectral ratio method along 3,852 ray paths between 708 aftershocks of the 2015 earthquake and a network of 22 broadband seismometers. We then used these measurements in a 3D nonnegative least-squares regression inversion to determine the subsurface QS structure. We identify high QS anomalies on the Nazca – South America subduction interface in the Illapel region which correlate spatially with areas of high frequency seismic radiation and high coseismic slip from the 2015 earthquake, and we interpret these anomalies as asperities. We also observe elongated bands of alternating high and low QS on the plate interface which match the orientation of normal fault fabrics on the surface of the subducting Nazca plate. We interpret the high QS bands as subducted horsts, which are materially strong and prominent, and behave as asperities; low QS bands are subducted grabens full of more attenuating subducted sediments and underplated continental material. These structures guided the rupture of the Illapel earthquake: nucleation occurred in a region of moderate QS (∼400) within a graben, but propagated to two adjacent high QS (> 1000) asperities which were critically stressed to near-failure, one on a horst directly up dip and another on a horst directly down dip. Both of these asperities subsequently failed, causing a cascading rupture which characterized the great earthquake.

SummaryOver the past few decades, the global number of dams has increased substantially. Water impounded behind these dams, resulting in elevated crustal pore pressure and altered stress distribution around reservoirs, could potentially trigger or suppress the failure of nearby faults, leading to transient changes in seismicity. In this study, we analyze 14 years (2006-2019) of spatiotemporal seismicity recorded by a dense local network in the Gotvand area (SW Iran), covering about 5.5 years before and 8.5 years after impoundment. The initial catalog, comprising over 48,000 relocated earthquakes, was reduced to 6,464 background events by adopting a 3D ETAS declustering model with a cutoff magnitude of 1.3. We analyze the spatiotemporal background seismicity pattern in the Gotvand area in comparison with calculated reservoir-related spatiotemporal stress changes relative to the initial stress state before water impoundment, approximating the Gotvand reservoir by 726-point sources covering the reservoir surface. We find that following the initiation of impoundment, the local background seismic activity slightly increased during the impounding, pointing to induced/triggered seismicity. However, most importantly, the impoundment of Gotvand lake has altered the spatial seismicity patterns, leading to a notable reduction in seismic activity in the central area of the reservoir, which is in agreement with the calculated negative Coulomb stress changes in the same area. Using the Coulomb-Rate-and-State seismicity model, we find that the spatiotemporal seismicity response due to the calculated stress changes is consistent with the observations.

SummaryWe present a statistically and computationally efficient spectral-domain maximum-likelihood procedure to solve for the structure of Gaussian spatial random fields within the Matérn covariance hyperclass. For univariate, stationary, and isotropic fields, the three controlling parameters are the process variance, smoothness, and range. The debiased Whittle likelihood maximization explicitly treats discretization and edge effects for finite sampled regions in parameter estimation and uncertainty quantification. As even the ‘best’ parameter estimate may not be ‘good enough’, we provide a test for whether the model specification itself warrants rejection. Our results are practical and relevant for the study of a variety of geophysical fields, and for spatial interpolation, out-of-sample extension, kriging, machine learning, and feature detection of geological data. We present procedural details and high-level results on real-world examples.

SummaryThe Surface Water and Ocean Topography (SWOT) altimeter mission provides high-resolution sea surface heights (SSHs), which can be used to retrieve high-precision and high-resolution marine gravity fields. However, data acquisition from the SWOT satellite’s Ka-band Radar Interferometer (KaRIn) is characterized by a 20-km nadir gap and a seam gap. Interpolation methods were employed to fill the nadir and seam gaps in the SWOT/KaRIn SSHs. The 20-km nadir gaps were interpolated using SWOT/KaRIn SSHs from flanking swaths on both sides. To reduce temporal decorrelation, seam gap interpolation was performed separately for each pass using SWOT/KaRIn SSHs from adjacent passes. The vertical gradient of gravity anomaly (VGGA) model was inverted from the interpolated SWOT/KaRIn SSHs using least-squares collocation and the remove-restore methods. The study region was selected as the Philippine Sea, with SWOT/KaRIn SSHs from cycle-02 serving as the experimental dataset. SWOT/KaRIn SSHs were processed using four schemes: Kriging interpolation, Cubic Spline interpolation (CSI), mean sea surface (MSS) interpolation, and no interpolation. The VGGA models inverted from the respective processed SWOT/KaRIn SSHs are hereafter referred to as Kriging_SWOT_VGGA, CSI_SWOT_VGGA, MSS_SWOT_VGGA, and NO_SWOT_VGGA. The results show that the NO_SWOT_VGGA model exhibits distinct strip artifacts, whereas the interpolated VGGA models eliminate these artifacts. All interpolated VGGA models exhibited superior consistency with the reference model compared to the NO_SWOT_VGGA model. Among them, the Kriging_SWOT_VGGA model achieved optimal consistency. Therefore, this study on interpolation methods for SWOT/KaRIn SSHs provides a novel methodological framework for the recovery of high-precision marine gravity fields from SWOT observations.

SummaryWe have analyzed direct S-waves of teleseismic earthquakes to investigate seismic anisotropy parameters, i.e., fast polarisation direction (FPD or φ) and splitting time delay (STD or δt) beneath southwestern Tibet (around the Karakoram Fault), that enable us to comprehend the upper mantle dynamics of the study region. To achieve this aim, we employ the Reference Station Technique, which is proven to be insensitive to source-side anisotropy; hence, it permits the use of teleseismic direct S-wave signals in splitting measurements. A total of 1,624 high-quality direct S-wave splitting measurements were obtained from 145 earthquakes (M ≥5.5) within an epicentral distance of 30○ to 90○, recorded at 31 temporarily deployed seismic stations of the Y2 network. We have found STDs ranging from 1.1 s to 1.8 s, indicating a significantly anisotropic upper mantle underneath the study region. Our splitting measurements reveal predominantly E-W FPDs in the western part of the study region, with a slight shift to the ENE-WSW direction in the eastern section. A comparison of our direct S-wave derived splitting measurements with prior SKS splitting measurements indicates a largely analogous pattern at most seismic stations. The seismic stations (WT04, WT05, WT11, and WT18), which previously lacked SKS-derived seismic anisotropy, are now complemented with new measurements with clear anisotropic signatures. Nearly E-W oriented FPDs that exhibit an oblique variation to the main strike of the southeastern segment of the Karakoram Fault (KKF) can be explained by the eastward movement of upper mantle material beneath southwestern Tibet. The significant discrepancies between the orientation of FPDs and the strike direction of KKF imply that the fault is not a lithospheric-scale fault but rather is confined to crustal depths. Integrating surface deformation derived from geodetic measurements (e.g., global positioning system data) and plate motion vectors of the Indian and Eurasian plates with splitting parameters indicates that the deduced deformation patterns result from both lithospheric deformation and sub-lithospheric mantle dynamics. The FPDs exhibit a significant deviation from GPS data, signifying a decoupling of crust and upper mantle materials beneath the study area. This suggests that mantle deformation in southwestern Tibet operates in a manner that is distinct from that of crustal deformation. Finally, our novel splitting measurements, enhanced by a greater number of direct S-wave data, provide new insights into the deformation of the upper mantle in the region, elucidating the mechanisms that have shaped the plateau over geologic millennia.

SummaryNormal modes or whole Earth oscillations serve as vital tools in the investigation of the Earth and other planets, providing large-scale seismic constrains on their globally averaged interior structure. Recently, Lognonné et al., 2023b reported their detection on Mars for the largest recorded event S1222a. However, the low signal-to-noise ratio in the normal mode frequency domain makes detection of normal modes very difficult, necessitating comprehensive data processing to extract the normal mode frequencies including deglitching, phase cross correlation, multi-tapering and phasor-walkout techniques. Here, we show that normal mode spectra for event S1222a depend on the details of deglitching technique, producing different results using either instrument response deconvolution (MPS) (Scholz et al., 2020) or the novel time-frequency polarization filtering (Brinkman et al., 2023). Furthermore, the plethora of Martian seismic models published since the InSight mission also generate strongly varying synthetic spectra, which complicates the identification of individual modes. Here, we developed a different approach in which we use spectra for the S1222a event and data for the seismically ‘quiet’ day prior to the event. We compare these observed spectra to synthetic spectra for existing seismic models for Mars in order to verify the detection of normal modes and identify models which best fit the spectra. Our research incorporates a range of different filtering and deglitching sequences and eleven post-InSight Martian seismic model families to assess their effects on both observed and synthetic spectra. We find that, the synthetic spectra consistently yield more overlapping peaks with the S1222a data spectra, than with spectra from the quiet day before S1222a. Through overlap analysis, we identify a preferred crustal structure, with an average Vs of 3.25 km/s, that aligns well with current geophysical and seismic estimates of the global Martian crust. We also find that the detection of normal modes is improved by stacking the data of S1222a with the two additional events S1000a and S1094b. On the other hand, stacking of all available long period seismic data of InSight to detect the continuous activation of normal modes by Martian atmosphere proved less effective. Nevertheless, we find consistent spectral peaks across the various data sets and stacking methods. These findings indicate that normal modes are detectable on Mars. It highlights the need for future deployment of seismometers on Mars, the moon and other planets to continue the hunt for normal modes and improve our understanding of the internal structure of terrestrial planets and moons.

SummaryAzimuthal anisotropy of the upper mantle, resolved from seismic records, sheds light on the characteristics of lithospheric deformation and asthenosphere flow. However, our understanding of the three-dimensional structures of P-wave azimuthal anisotropy in the upper mantle remains limited, due to the lack of accurate and high-resolution tomography methods. Motivated by this necessity, we present a novel full-wave inversion method to simultaneously retrieve the three-dimensional structures of P-wave velocity and azimuthal anisotropy of the upper mantle. The inversion is parameterized by imposing five elastic coefficients representing Pn-wave azimuthal anisotropy on the isotropic elastic tensor based on the weak anisotropy assumption. We verify the accuracy of the full-wave sensitivity kernels of Pn waveform cross-correlation delay times to the azimuthal anisotropic coefficients. Synthetic inversions demonstrate that our method effectively resolves the anomalies of both P-wave velocity and azimuthal anisotropy including 2ψ and 4ψ terms. Moreover, we analyze the model resolution and tradeoffs using the Hessian-vector product efficiently computed by the Scattering-Integral method. This study provides a novel and powerful physics-based tool to reveal accurate and high-resolution three-dimensional structures of azimuthal anisotropy in the upper mantle, which will facilitate our understanding of the modes of lithospheric deformation and mantle convection.

SummaryFull waveform inversion (FWI) updates the velocity model by minimizing the discrepancy between observed and simulated data. However, incomplete seismic acquisition can introduce errors that propagate through the adjoint operator, affecting the accuracy of the velocity gradient and reducing the convergence accuracy and speed. To mitigate the influence of acquisition-related noise on the gradient, we employ a convolutional neural network (CNN) to extend the velocity representation and refine the velocity model before forward simulation, thus reducing gradient noise and providing a more accurate velocity update direction. The same data misfit loss is used to update both the velocity and the network parameters, forming a self-supervised learning procedure. Here, the CNN acts as a dynamic velocity conditioner that is optimized to help fit the data. In this method, the velocity representation is extended (VRE) by combining a neural network with conventional grid-based velocities. Thus, we refer to this general approach as VRE-FWI. Synthetic and real data tests demonstrate that the proposed VRE-FWI achieves higher velocity inversion accuracy compared to traditional FWI, with only a marginal additional computational cost ∼1%.

AbstractAnisotropy of magnetic susceptibility (AMS) analysis is widely used as an efficient petro-fabric tool to infer magma flow patterns within dikes. However, interpretations of magnetic fabric often get complicated by the occurrence of anomalous (intermediate/inverse) fabrics oriented normal to the dike plane, which may lead to uncertainty, unlike the straightforward normal fabrics along the intrusion plane. In this article, we present a detailed rock-magnetic and magnetic fabric study of India’s ~116 Ma old Salma dike, which is the most prominent and longest dike related to the early Cretaceous Rajmahal Trap (RT) volcanism. A joint analysis of in-phase and out-of-phase anisotropy of magnetic susceptibility (i.e. ipAMS and opAMS) and anhysteretic remanent magnetization (AARM) fabrics allowed us to identify the different sources of the observed fabric and subfabrics. Rock-magnetic analyses suggest that the magnetic mineralogy consists of at least two types of titanomagnetite with varying Ti-content. FORC suggests that the dike is dominated by PSD grains with varying influence from SD grains. The ipAMS fabric is primarily carried by SD grains of low-Ti titanomagnetite, while opAMS is governed by larger MD/PSD titanomagnetite with higher-Ti content. Results indicate that anomalous, intermediate-type ipAMS fabrics, particularly along the dike margins, are caused as a combined effect of SD grains, late-stage crystallization, and mild high-temperature oxidation. At the dike center, post-emplacement alteration, intense exsolution of the primary titanomagnetite, and magma backflow are responsible for anomalous fabrics. In contrast, the majority of normal ipAMS, opAMS, and AARM fabrics are coaxial, providing a reliable record of magma flow as confirmed from the long axes trend distribution of the plagioclase laths along the dike plane. These fabrics reveal a dominant subvertical magma flow direction during emplacement, indicating magma ascent from depth. Even though the magnetic fabrics do not unequivocally constrain the deeper processes underneath Moho, they are compatible with the idea of a subcrustal magmatic layer beneath the region, potentially originating from decompression melting or the influence of the Kerguelen plume, to be the feeder source.

SummaryLocated at the easternmost passive margin of the Eurasian Plate, the southeastern Korean Peninsula shows geological signatures consistent with Miocene backarc opening associated with the Pacific Plate subduction. The region comprises two contrasting crustal blocks—the Early Cretaceous Gyeongsang Basin and Miocene Yeonil Basin (YB)—and hosts multiple fault systems that record both extensional and contractional deformation, providing an ideal setting to investigate crustal evolution along a passive margin. Motivated by this complex setting, we performed high-resolution P-receiver function imaging using a dense broadband seismic network. Our results reveal two Moho offsets: a western offset from 33 to 28 km, and an eastern offset from 28 to 26 km, coinciding with major fault zones and likely reflecting localized crustal thinning and subsequent reactivation. Crustal anisotropy, inferred from changes in fast-axis orientations, varies spatially, with Miocene fossil anisotropy in the GB and both fossil and present-day stress-induced anisotropy in the YB. Variations in P-to-S velocity ratio (VP/VS) reflect compositional heterogeneity and fault-related fracturing. Large earthquakes (M ≥ 4) occurred in low-VP/VS zones associated with relatively rigid and possibly locked crustal segments, while high-VP/VS regions coincide with zones of crustal weakening and microseismicity. Our findings suggest that extension-related deformation and inherited structural heterogeneity are preserved within the crust of this fossil backarc system, linking past tectonic processes to present-day structure and seismicity.

SummaryThe quantification of frozen and unfrozen water content in porous media is essential for understanding hydrological, thermal, and mechanical processes in cold regions. Petrophysical joint inversion (PJI) frameworks that integrate seismic and electrical data offer promising tools for resolving water and ice distributions but are often limited by simplified petrophysical models. Here, we extend a PJI framework by incorporating a temperature-dependent relationship that accounts for both electrolytic and surface conduction. Using synthetic data, we show that this formulation improves modelling of frequency-dependent resistivity and enhances estimates of water content and interfacial conductivity quantified by cation exchange capacity. Our results highlight the critical role of temperature in controlling subsurface electrical properties and demonstrate that neglecting these effects can lead to substantial errors in the ice and water estimates. The extended PJI framework provides a physically consistent basis for geophysical imaging of water phase dynamics in partially frozen systems, with broad applicability to cold-region hydrology and seasonally or perennially frozen environments.

SummaryWe simulate an instantaneous time mirror (ITM), i.e., a rapid short duration change in elastic material properties, using numerical experiments in time-varying isotropic elastic media. Our implementation in the seismic wave propagation software SeisSol is based on high-order discontinuous Galerkin discretization with ADER time stepping. We develop an eigenvector-based analytical solution for time interfaces for general linear hyperbolic wave systems and apply it to analyze the energy balance at time boundaries and ITMs. The energy increases for all intermittent medium changes for all impedance scaling factors. Our numerical implementation is validated against these analytical solutions and achieves high-order convergence. Its accuracy is further corroborated by estimates of reflection and transmission coefficients and observed frequency shifts across time boundaries, and by acoustic wave speed estimates obtained from focal spots associated with ITM-generated converging P waves that are consistent with theoretical predictions and ground truth values, respectively. We use the ITM implementation to simulate the partitioning of seismic body waves excited by a point source in a spatially homogeneous elastic full space. The response to an intermittent short change in the elastic parameters yields a diverging and converging P and S wavefield. A systematic scaling of the elastic parameters is then used to steer independent ITM reflections of either P or S waves. Numerical ITM solutions as developed here can be used to synthesize converging wavefields in seismic imaging applications, and more generally to analyze the behavior and manipulation of seismic wavefields in space-time varying media.

SummaryAccurate characterization of the magnitude-frequency distribution of seismicity, and its associated uncertainties, is essential for seismic hazard assessment. This distribution is commonly described by the Gutenberg–Richter (GR) relation, parameterized by the b-value, which has been identified as a potential proxy for investigating many spatiotemporally varying Earth phenomena. Estimating the spatiotemporal variability of b-values often requires windowing, forcing a trade-off between resolution and statistical reliability. New probabilistic methods circumvent this by inferring both the number and locations of change points directly from earthquake catalogs. Nevertheless, accurately determining the b-value remains difficult because the GR relation only holds over a limited range of magnitudes. This research develops a general statistical model to address several methodological challenges in estimating the magnitude-frequency distribution of observed seismicity, including variations in space or time. The approach simultaneously solves for the b-value and magnitude-range limits. This avoids potential bias due to inaccurate manual truncation of earthquake catalogs. The model considers the entire observed catalog and parameterizes the decay of the distribution at both low and high magnitudes. Consequently, robust uncertainties in estimated b-values reflect uncertainty in the range of magnitudes over which the GR relation is observed to be valid. Importantly, spatiotemporal variations in the parameters that define the magnitude range are considered to be independent from the b-value, as we assume the physical factors that influence the GR relation are independent of the factors that limit the observed earthquake catalog. We demonstrate this methodology through application to simulated and observed earthquake catalogs. In particular, the value of our approach is highlighted through application to observed records of induced seismicity associated with fluid-injection operations in western Canada. Our results demonstrate accurate b-value estimates and associated uncertainties. Furthermore, the additional parameters that define the magnitude range serve as proxies for other factors including seismic network performance, recording duration, potential geometric limitations on earthquake size, and potential injection characteristics (in induced seismicity cases). Our approach also allows for the investigation of how these other factors may vary in space/time. Results from this work contribute to rigorous propagation of accurate b-value estimates, including uncertainties, into subsequent analyses such as seismic hazard models and regulatory protocols that are applied to industrial activity.