Geophysical Journal International

AbstractSeismic inversion translates seismic data into subsurface elastic property models, enabling geophysicists to better understand underground rocks and fluids. Due to the inherently ill-posed nature of this inverse problem, accurately capturing the uncertainty associated with the solution is essential for reliable interpretations. Traditional Bayesian inversion methods, such as Markov Chain Monte Carlo (MCMC) and Laplace approximations, have been employed for this purpose but face significant limitations in terms of scalability and computational efficiency for large-scale problems. Combined with deep learning, Variational Inference (VI) has emerged as a promising alternative, striking a balance between computational efficiency and flexibility (i.e. the ability to approximate complex posterior distributions). However, selecting an appropriate proposal distribution remains a key challenge, as it directly influences the quality of the estimated posterior distribution. In this study, we extend IntraSeismic, an implicit neural representation (INR)-based framework for seismic inversion applications, to Bayesian inversion using VI with different parameterizations of the proposal distribution. We introduce two methods: B-IntraSeismic (BIS), which uses a mean-field Gaussian proposal, and B-IntraSeismic with Conditional Normalizing Flows (BIS-Flow), which utilizes a mean-field unparameterized proposal distribution to better capture deviations from Gaussianity in the posterior distribution. These methods are evaluated on a synthetic dataset (Marmousi) and two field data (Volve and Sleipner). Our results indicate that both BIS and BIS-Flow can accurately capture structural details and produce high-resolution mean models and standard deviation maps. BIS-Flow is also shown to be able to model complex posterior distributions, offering a more comprehensive characterization of uncertainty while maintaining computational feasibility.

AbstractThe phenomenon of elastic wave conversions, where acoustic, pressure (P-) waves are converted to elastic, shear (S-) waves, and vice-versa, is commonly disregarded in seismic imaging. This can lead to lower-quality images in regions with strong contrasts in elastic parameters. While a number of methods exist that do take wave conversions into account, they either deal with P- and S-waves separately, or are prohibitively computationally expensive, as is the case for elastic Full-Waveform Inversion. In this paper an alternative approach to taking converted waves into account is presented by extending Full Wavefield Migration (FWM) to account for wave conversions. FWM is a full-wavefield inversion method based on explicit, convolutional, one-way propagation and reflection operators in the space-frequency domain. By applying these operators recursively, multi-scattering data can be modelled. Using these operators, the FWM algorithm aims to reconstruct the reflection properties of the subsurface (i.e. the ‘image’). In this paper, the FWM method is extended by accounting for wave conversions due to angle-dependent reflections and transmissions using an extended version of Shuey’s approximation. The resulting algorithm is tested on two synthetic models to give a proof of concept. The results of these tests show that the proposed extension can model wave conversions accurately and yields better inversion results than applying conventional, acoustic FWM.

SummaryGlobal seismology mainly uses seismic waves propagating in the sagittal plane along the great circle path (GCP). However, heterogeneities in the mantle laterally deviate the path of seismic signals, which arrive out-of-plane (OOP) at arrays of sensors at teleseismic distances. Detection and back-projection of these signals have, in the past, provided independent evidence for the location of distant subducted slabs in the deep mantle, complementing global tomographic imaging. To infer physical properties of these subducted slabs, 3D waveform modeling of OOP waves for a finite-thickness slab is needed but still missing. In this study, we conduct a series of synthetic tests using a spectral element solver. We test the detectability of OOP signals and, by progressively adding complexities, we evaluate to which extent these signals can be used to infer physical properties of the modeled slab. We carry out three-component array analysis and investigate focal mechanism dependency. Our results show that the transverse component might be the best candidate for such studies, also for P-to-P OOP signals. Vertical and radial component recordings are usually dominated by P-SV energy arriving from the earthquakes along the GCP, which masks possible OOP signals. Contrary, the transverse component filters out any P-SV energy arriving directly from the source and, owing to its intrinsic directionality, allow for higher resolution measurement of P-to-P OOP signals. This is especially the case prior to the arrival of the S-wavefield. We pick a series of OOP arrivals which are back-projected using a multi-phase trial-and-error approach, that is considered successful only when different OOP seismic phases converge to the modeled (true) structure. We retrieve the location of the slab, its bottom and top edges, and its thickness in the lower mantle. These inferences are tested against varying topography, orientation and size of the modeled slab. The insights gained with modeling are confirmed with real data examples, supporting higher resolution mapping of 3D mantle structure based on OOP seismology.

SummaryWe investigate the lithospheric structure of the Southwest Iberian margin along an active seismic profile southwest of São Vicente Cape, ranging from the southern Tagus Abyssal Plain to the westernmost part of the Gulf of Cadiz. This profile, approximately 320 km long, intersects almost perpendicularly three major thrust faults: the Tagus Abyssal Plain, Marquês de Pombal and Horseshoe faults. The crustal structure, derived from spatially coincident wide-angle seismic (WAS) and multichannel seismic (MCS) data, was validated and constrained using gravimetric data. Joint travel-time inversion of refracted phases identified in WAS and reflected seismic phases from both WAS and MCS records were used to build a detailed two-dimensional P-wave velocity (Vp) structure. The resulting model reveals a Vp distribution with abrupt lateral velocity and structural variations, characterized by a rugged basement top and sharp changes in crustal thickness. Three main lithospheric domains consisting of continental, oceanic, and exhumed mantle affinity were identified from south to north. The travel-time inversion of the deepest reflected seismic phases reveals four major southeast-dipping reflectors, likely corresponding to major regional thrust faults with significant seismic and tsunamigenic potential. Integrating the modelled and interpreted seismic results with the locations of recent well-constrained earthquakes suggests that the Marquês de Pombal and Tagus Abyssal Plain extend deeper than previously thought, with fairly high seismic activity in the deep levels. This has significant implications for their seismogenic potential and should be taken into account for accurate assessment of seismic hazards in the region.

SummaryFirst-arrival traveltime tomography (FATT) is a widely used method for characterizing near-surface velocity structures in geotechnical engineering and resource exploration. We introduce an improved version of FATT, named first-arrival adjoint tomography (FAAT), which involves the joint inversion of first-arrival absolute traveltimes and differential traveltimes. Unlike absolute traveltimes, differential traveltimes, derived from common sources or receivers, offer heightened sensitivity to fine-scale structures near the receivers or sources, respectively. This dual sensitivity makes FAAT particularly effective in imaging highly heterogeneous media. However, the proximity of rays associated with differential traveltimes can lead to the instability of the inversion. To overcome this challenge, we simultaneously incorporate absolute traveltimes and differential traveltimes to update the velocity model. This approach improves the stability of the inversion process, leading to improved resolution of inverted results. We specifically employ the fast-sweeping method to solve the factored eikonal equation, providing robust solutions in models with complex geological structures. Furthermore, we address the inverse problem by computing the gradient of data misfit using the efficient adjoint-state method. Through numerical testing, we validate the effectiveness of FAAT in comparison to that using only absolute or differential traveltimes. Finally, we apply the proposed FAAT method to near-surface site characterization at Bishan-AMK Park in Singapore. Compared with FATT and validated against borehole data, FAAT demonstrates its ability to establish more reliable velocity models, revealing finer details and substantially improving geological interpretation.

SUMMARYNumerical models are a powerful tool for understanding and predicting the impact of landslide-generated tsunamis. We use here the HySEA code which incorporates a multilayer structure and non-hydrostatic pressure to simulate the tsunami generated by a potential submarine landslide located offshore of the Mayotte island. The island is surrounded by a lagoon and steep slopes weakened by the 2018 seismo-volcanic crisis. The influence of the input parameters and of the model assumptions is shown to change by a factor 2 the predicted maximum water free surface elevation, velocity and maximum inundation depth. This demonstrates the need of using numerical models for building local scaling laws to relate tsunami and landslide properties. Our results highlight the necessity of incorporating high-resolution bathymetry, in-depth variations through multilayer modeling and relevant landslide rheology to accurately predict tsunami impact. In case of strong topography variations as in Mayotte, using 4 layers seems to be a good compromise between accuracy and computational cost. Accounting for these effects would enable to refine hazard maps by identifying safe and high-risk coastal zones and to improve wave arrival time estimates, thus reducing tsunami-related risks in regions like Mayotte.

SummaryUnder compression in the brittle regime, rocks fail as fracture interaction, propagation and coalescence produce a process zone of high strains that forms the nucleus of the fault that eventually traverses the rock core. Early analyses suggested that this process zone extends through relatively intact rock close to the peak stress, and that the deformation that develops earlier in loading does not guide the location of this process zone. Here, we assess the predictability of the location of the process zone in triaxial compression experiments on Darley Dale sandstone. We develop supervised machine learning models to predict whether a particular location in the rock core will eventually host high shear and dilative strain, and thus eventually reside within the process zone. We calculate the local incremental strain tensors throughout the rock cores using digital volume correlation of X-ray tomograms captured during the experiments. We find that the models produce accuracies of 0.73 on average, indicating that it is possible to predict the location of the process zone using observations of the strain field preceding macroscopic failure. The average accuracy increases to 0.84 when the data are restricted to later in the experiment because the strain field becomes more similar to the final strain field that hosts the process zone. The models primarily rely on the shear strain field to predict the location of the process zone, likely because it is more spatially persistent than the dilative strain. The local neighborhood of the strain field is helpful for predicting whether a location will eventually host the process zone only when the data are restricted to one experiment. When the data includes several experiments, the models primarily rely on global statistics of the strain fields. The varying correlation lengths of the dilative and shear strain field of different experiments help explain this result.

AbstractActive tectonic movements and surface deformation are observed in the eastern Tibetan Plateau. Understanding variations in crustal thickness and the deep Moho interface is crucial for elucidating the expansion of the Tibetan Plateau. This study utilizes InSAR to derive vertical surface deformation and applies loading corrections based on Green's function method. Additionally, satellite gravity data are used to separate hydrological and tectonic signals to infer changes in the Moho interface. Our results indicate that the regionally averaged loading effects, estimated using localized Green's functions (LGF), account for approximately 16.5 per cent of the InSAR-derived vertical displacement field. This contribution exhibits significant spatial variability, exceeding 100 per cent in regions with strong hydrological activity. The loading calculation is highly sensitive to Earth model: the relative difference between the load displacement obtained using the local Green's function and that obtained with an average Green's function reaches 48.2 per cent. After applying loading corrections, a more accurate Moho uplift rate of –8.2 ± 3.1 mm/a is obtained. The findings support the conclusion that the Moho interface rises in the southern region of the study area, with crustal thinning, while the Moho surface sinks in the northern region, with crustal thickening.

AbstractIn recent years, offshore tsunami observation networks equipped with ocean bottom pressure gauges (OBPGs), such as S-net, DONET, and N-net, have been deployed around Japan, enabling real-time collection of high-quality tsunami data near the source. These networks make it possible to estimate the spatiotemporal variation of the tsunami wavefield using a data assimilation approach, and to predict coastal tsunamis from the initial or current tsunami wavefield. This study proposes a novel tsunami data assimilation method that uses physics-informed neural networks (PINNs) to estimate tsunami wavefields from the observed OBPG data. The neural network was optimised by minimising the sum of the data loss, which quantifies discrepancies from the tsunami data, and the physical loss, which quantifies the satisfaction of the linear long wave equation. This was performed to ensure that the estimated results are consistent with both the observed data and the physics of tsunami propagation, even when there are limited observational data and significant noise. We first validated the effectiveness of the proposed method using synthetic S-net OBPG data from the 2011 Tohoku-oki earthquake (Mw 9.0) tsunami. The results confirmed that by using both data and physical constraints in the PINN optimisation, the PINN could adequately assimilate the spatiotemporal distribution of the tsunami wavefield from OBPG data, even for predictions outside the network coverage area. The predicted tsunami waveforms at the coastal stations, computed from the estimated initial wavefield, showed good agreement with the actual waveforms. Next, we conducted an experiment using actual S-net OBPG data from the 2016 Fukushima-oki earthquake (Mw 6.9) tsunami. The initial tsunami source estimated by PINN was in good agreement with other studies based on waveform inversion, although the maximum source amplitude and maximum coastal tsunami heights were underestimated. We also conducted an experiment using N-net OBPG data from the 2024 Hyuganada earthquake (Mw 7.0) tsunami. The PINN could accurately estimate the initial tsunami source, even though the tsunami source of this event was located outside the N-net coverage area. Finally, we have shown that incorporating tsunami observations over time into the iterative optimisation of the PINN model allows for accurate and efficient tsunami data assimilation.

AbstractUnderstanding processes in the Critical Zone requires reliable information about the vadose-zone aquifer, its geometry, and spatial variability. Commonly, such information is obtained from boreholes, yet large areas might render their application prohibitively expensive. Additionally, limited geological a-priori information might bias the interpretation due to lateral geological changes smaller than the borehole sampling scale. The transient electromagnetic method (TEM) has emerged in the last decades as a well-suited method to efficiently investigate the subsurface, as required for many hydrogeological applications. The interpretation of TEM measurements relies mainly on deterministic inversions, offering only a limited insight on the uncertainty of the subsurface model. Uncertainty quantification, however, is essential for integrating TEM results into hydrogeological models. Hence, we propose a combined approach using both deterministic and stochastic inversion of TEM soundings to investigate the uncertainty of shallow (< 40 m) aquifers. Current stochastic approaches for TEM data rely on Markov chain Monte Carlo algorithms, which have to be run from scratch for each individual sounding. Alternatively, machine learning approaches, such as Bayesian Evidential Learning (BEL), can be much faster because they do not require retraining for every new data set. Hence, we investigate, in particular, the application of a single, common prior model space instead of multiple, individual prior model spaces to directly estimate the uncertainty of multiple TEM soundings. To this end, we combine forward modelling routines with the stochastic inversion approach BEL1D and assess our approach using both field data and numerical experiments.

SUMMARYThe Fennoscandian earthquake catalogue (FENCAT) assembles data on the natural seismicity in Fennoscandia, Northern Europe. We present an updated and standardized version of the catalogue originally published in the early 1990s. New instrumental data are recorded by the seismic networks of Denmark, Estonia, Finland, Norway, and Sweden, and analyzed by the Geological Survey of Denmark and Greenland, the Geological Survey of Estonia, the University of Helsinki in Finland, the University of Bergen and the NORSAR research foundation in Norway, and Uppsala University in Sweden. The updated catalogue provides the available earthquake parameters in a brief, user-friendly version: origin time, source coordinates, focal depth, macroseismic data (maximum intensity and radius of the area of perceptibility), up to three observed magnitudes, seismic moment estimate, and a standardized moment-related magnitude, mW(HEL), for each event. The standardized magnitude is defined in this paper and its relation to other magnitude scales is provided. Suspected non-earthquakes (e.g. frost events, explosions, human-induced events) have been removed. The standardized event magnitudes range from mW(HEL) -1.0 to 6.2. To enable the usage of earthquake data in a large variety of seismological, geological and earthquake engineering investigations, the data are not truncated at the low magnitude end.The updated catalogue, FENCAT (2021), contains about 23 000 earthquakes for the period 1467–2021 in an area bounded by 54–75°N latitudes and 0–45°E longitudes. The completeness and quality of the earthquake solutions is best within the areal coverage of the above-mentioned networks.

AbstractWhile modern thermal convection in rocky planets is controlled by a slow solid-state creep flow, the earliest stages of terrestrial planets likely experienced turbulent flow during which their silicate envelope was fully molten, usually called magma ocean. The main parameter separating the two regimes is the Prandtl number (Pr), which is so high for mantle convection to be usually assumed infinite, whereas magma oceans are characterized by Pr on the order of 1. We compared the results of isoviscous convection simulations performed with three codes: (GAIA, TLBM, StreamV). These codes are based on different numerical formulations and were used for modeling convection with Pr ranging from 1 to 1000, while exploring different convection intensity by varying the Rayleigh number (Ra) from 104 to 106. GAIA (Generic Automaton for planetary Interior Analysis) is a Finite Volume fluid flow and energy solver for the Navier-Stokes equations across arbitrary geometries. TLBM (Thermal Lattice Boltzmann Method) solves the mesocale momentum and energy distribution densities for colliding particles on a discrete lattice. StreamV is a Eulerian-Lagrangian Finite Volume code that solves the Navier-Stokes equations under the Boussinesq approximation. The codes are compared over 24 different simulation setups, analogue to the classical Blankenbach infinite Pr benchmark (Blankenbach et al., 1989a), but extending it to finite Pr and to two types of boundary conditions, free-slip and no-slip. We show that the results of the three codes are generally in good agreement, and discuss differences. Finite Pr solutions show a much richer dynamics varying from stable steady-state solutions, to oscillatory and chaotic ones, and converging to infinite Prandtl number solution for increasing values of Pr for larger Ra: Pr ≥ 100 is sufficient for Ra = 105 but Pr ≥ 1000 is required for Ra = 106. Our results offer a robust set of solutions useful for testing future finite Prandtl number convection codes.

SUMMARYOur viewpoint highlights the challenges faced by women in the Induced Polarization (IP) community (and elsewhere), particularly the persistent gender imbalance in scientific workshops. This underrepresentation in leadership roles, presentations, and discussions reflects broader systemic biases in academia. By sharing the experience of the 7th IP workshop, where the organizing team made deliberate efforts to prioritize gender and diversity in organizing committees, recruitment, and session formats, we demonstrate how intentional actions can create a more inclusive, gender-balanced environment. This approach is crucial not only for the IP community but for all research communities. Emphasizing diversity and inclusion fosters a welcoming atmosphere that encourages participation from all individuals, regardless their background. In turn, this enriches the research process by incorporating diverse perspectives, driving innovation, and improving scientific outcomes. We aim to inspire other academic communities to actively promote diversity and inclusivity, recognizing that such efforts not only benefit underrepresented groups but elevate the success and relevance of science as a whole.

SummaryWe present a novel model-based seismic data redatuming method based upon time-domain data-assimilated wavefield (DAW) reconstruction. Seismic redatuming refers to the positioning of sources and/or receivers from the acquisition surface to a virtual datum level along with the computation of the dataset that would have been recorded by the virtual acquisition. Physically, DAWs approximate the true wavefields by least-squares time-reversal extrapolation of the recorded data in a background subsurface model. To this end, DAWs are first formulated as the sum of the background wavefields and an approximation of the scattered wavefields by the unknown contrasts between the true model and the background model. We estimate the scattered fields via the estimation of their volume scattering sources by fitting in a least-squares sense the differences between the recorded data and the simulated data in the background model, namely the restriction at receivers of the scattered fields. These underdetermined linear source problems involve two main steps: a multidimensional deconvolution of the scattered data by the data-domain Hessian, followed by the propagation of the deconvolved scattered data backward in time. Once the source problems have been solved, DAWs follow by solving the wave equation in the background model using the sum of the experimental sources and the scattering sources as extended sources. Finally, the redatumed datasets, whose accuracy depends on the precision of the background medium between the acquisition surface and the datum, are readily obtained by sampling the DAWs at arbitrary datum levels. The primary challenge in computing accurate DAWs lies in the multidimensional deconvolution of the scattered data, which requires solving a data-domain normal equation with preconditioned Krylov-subspace iterative methods, where each iteration requires one forward and one backward simulation. We demonstrate the accuracy of the method and discuss its computational cost with several benchmarks representing various experimental environments (onshore with weathered layers, shallow and deep offshore).

SummaryDistributed acoustic sensing (DAS) has emerged as a potential solution to the sparse instrumentation issue in the world’s oceans. DAS involves repurposing fibre optic cables into dense receivers. The spatial undersampling limits our understanding of fundamental oceanic processes, like ocean dynamics. We use long-term DAS recordings from Svalbard, Norway, over two roughly perpendicular fibre segments to analyse ocean surface gravity wave (OSGW) signals and gain additional insight into their dynamics. This fibre layout allows estimation of the angle of arrival for OSGW generated under different weather conditions, while the long-term recording allows one to study seasonal variations. By investigating different wind directions, we observe two sets of OSGW arrivals: swells generated by distant storms and waves generated by the local winds. The swells consistently originate from the south-west, whereas the wind-forced OSGW follows, more or less, the local wind direction. Moreover, we conduct a detailed analysis of the recorded swell waves by computing their origin time, great circle propagation distance, group velocity, incidence angle, location, and interference pattern. This yields important data that can be used to characterise local and distant Atlantic storms. Only one receiver system has been employed, and more receivers are needed to validate the results obtained here and gain additional insight into the OSGW signals recorded on DAS systems.

SummaryInversion of geophysical data usually exhibits strong non-uniqueness, arising from sparse data coverage, limited number of measurements, inherent nonlinearity of governing physical laws, noise, and other factors. Methods based on Monte Carlo sampling are commonly used to explore the posterior model distributions, but these approaches are computationally demanding. Variational inference (VI) provides an alternative by transforming a high-dimensional sampling problem into an optimization problem, thereby significantly reducing the computational time. However, conventional VI methods, which typically use simple distribution families, like Gaussians, to approximate the posterior, may lack flexibility necessary to capture the complexity of the posterior distributions. Normalizing flows (NFs), a type of deep generative models, address this limitation by transforming a simple initial distribution into a highly complex target distribution through a sequence of invertible and differentiable transformations. In this study, we develop an NF-based VI method and apply it to electromagnetic (EM) data. This approach allows for explicit integration of prior knowledge and reference models into the inversion process. Both synthetic tests and field applications on EM data demonstrate that NF-based inversion effectively recovers the posterior model distribution in a more efficient manner, while providing excellent data fitting performance. Unlike many other machine learning algorithms, NFs do not require a training set, making it highly transferable across various inversion problems with minimal adjustments. The proposed NF-based method offers a more robust and computationally efficient solution to uncertainty quantification and shows great potential to be extended to solve 3-D geophysical Bayesian inversions, a major challenge that the geophysical community has faced for decades.

SummarySoutheast Asia, bordered by significant tectonic plates such as the Indo-Australian, Pacific, and Philippine Sea Plates, is distinguished by its frequent tectonic activity and complex geological structures, making it one of the most dynamically evolving regions worldwide. In this study, we introduce a novel 3D P-wave velocity model of the upper mantle and transition zone in Southeast Asia using regional seismic traveltime tomography based on first-arrival data from the International Seismological Center. We employ an adjoint-state tomography method with normal-vectors independence to accurately invert for 3D velocities using a 1D reference model. Synthetic tests confirm the reliability of our model in delineating features of the subduction zone and the surrounding region. Our inversion results highlight distinct subducted slabs within the subduction zone and a pervasive low-velocity zone beneath Sundaland, which may be associated with lithospheric thermal weakening. Additionally, a mushroom-shaped low-velocity anomaly attributed to the Hainan mantle plume is identified beneath Hainan Island. The low-velocity anomaly observed beneath the western part of the Java Sea may be attributed to the combined effects of Sunda-Java slab subduction, lower-mantle flow through the Sunda Strait, and the influence of the Hainan mantle plume. Notably, beneath the Andaman Sea, we observe an east-west elongation of the northern Sumatra slab, potentially linked to the clockwise rotational opening of the Andaman Sea. Additionally, three potential rifts are identified beneath the subducting Sumatra-Java slab: beneath the Toba Volcano, the Sunda Strait, and the eastern segment of Java Island. Extensive high-velocity anomalies beneath the Philippine Islands and the South China Sea suggest a double-sided subduction process involving the Proto-South China Sea slab.

SummaryInversion of a given geophysical dataset cannot be complete without assessing the resolution and uncertainties associated with the model obtained. However, model appraisal may still be a challenging task from both a theoretical and a computational point of view. To tackle the problems of model estimation and appraisal, we introduce the Subtractive Optimally Localized Averages (SOLA) method to the geophysical electromagnetic community, through the example of linear inversion of induced polarization (IP) data. SOLA is a variant of the Backus-Gilbert method: it is computationally more efficient but also allows one to specify directly the target local averages of the Earth’s properties to be estimated, including their uncertainties. SOLA offers great flexibility in the construction of averaging kernels, via the design of target kernels, and direct control over the propagation of data errors into the local-average estimates. With SOLA we obtain a collection of i) local averages of the ‘true’ Earth model, accompanied with their ii) averaging kernels and iii) uncertainties. We investigated the performance of SOLA for the 2–D tomographic inversion of a field IP data set. The obtained chargeability model compares well with previous studies, and, most importantly, its resolution (the spatial extent of the averaging kernels) and uncertainties can be interrogated. We conclude that SOLA is a promising approach for geophysical-electromagnetic linear(ised) tomographies. In the case of IP inversion, to construct chargeability models and evaluate their robustness.

SummaryFault systems have geometrically complex structures in nature, such as stepovers, bends, branches, and roughness. Many geological and geophysical studies have shown that the geometrical complexity of fault systems in nature decisively influences the initiation, arrest, and recurrence of seismic and aseismic events. However, a vast majority of models of slip dynamics are conducted on planar faults due to algorithmic limitations. We develop a 3D quasi-dynamic slip dynamics model to overcome this restriction. The calculation of the elastic response due to slip is a matrix-vector multiplication in boundary element method, which can be accelerated by using Hierarchical Matrices. The computational complexity is reduced from O(N2) to O(Nlog N), where N is the number of degrees of freedom used. We validate our code with a static crack analytical solution and the SEAS benchmark/validation exercise from Southern California Earthquake Center. We further employ this method on a realistic fault system with complex geometry that was reactivated during the 2023 Kahramanmaraş–Türkiye doublet earthquakes, generating slip sequences that closely match real observations.

SUMMARYThe spherical harmonic coefficient Level-2 products of the Gravity Recovery and Climate Experiment (GRACE) mission are affected by north-south stripe noise. Toward this end, we have developed a new filter named Variational Mode Decomposition spatial (VMDS) filter that transforms the Equivalent Water Height (EWH) map derived from GRACE level-2 product into a one-dimensional sequence, which is then filtered by using variational mode decomposition. This approach overcomes the limitations of the singular spectrum analysis spatial (SSAS) filter, which well performs in the medium-frequency band but omits the high-frequency NSS noise. We thus put the VMDS filter behind the SSAS filter to utilize the good performance of SSAS in the medium-frequency band and thus propose a combined filter termed SV. A closed-loop simulation demonstrates the better ability of SV to suppress NSS noise and preserve signal at the grid scale compared to the SSAS filter. In the real-world scenario, the SV solution achieves a noise level (46.68 mm of EWH) below that for SSAS and DDK7 solutions (53.52 and 53.68 mm of EWH, respectively) over the ocean at low latitudes. Moreover, the well-documented water level of Lake Victoria and the well-modeled coseismic gravity change of the Mw9.2 2004 Sumatra-Andaman earthquake demonstrate that the SV filter efficiently preserves localized mass evolutions while suppressing north-south stripe noise. Such short-wavelength signals usually miss in highly filtered spherical harmonics (e.g. DDK5 and DDK6) solutions or are significantly inconsistent for various mass concentration solutions.