Geophysical Journal International

SummaryWe present a probabilistic framework for evaluating earthquake forecasting models that use an alarm-based approach. In this approach, alarms are triggered by specific precursor signals. In a previous paper we compared such models and two ensemble models combining them in additive and multiplicative mode, with the ETAS (Epidemic Type Aftershock Sequence) forecasting model, which is defined in a probability-based approach, by making the latter to issue an alarm when the expected rate exceeds a predefined threshold. In this work we compare the alarm-based models with the ETAS and with another probability-based model, EEPAS (Every Earthquake a Precursor According to Scale) previously applied to Italy, using the testing procedures developed for probability-based models within the Collaboratory Study for Earthquake Predictability (CSEP) initiative. To do that, for the four alarm-based models, we compute empirical probabilities (frequencies) of Mw ≥ 5.0 earthquakes in Italy, inside and outside alarm time intervals issued by such models from 1990 to 2011. We then compare pseudo-prospectively the forecasting ability of all six models, by applying the CSEP tests on the time interval from 2012 to 2023. We found that the evaluation method used has a strong impact on the ranking of model performance. Probabilistic models like ETAS and EEPAS tend to score better under the CSEP testing framework whereas alarm-based models generally outperform probability-based ones when assessed using alarm-based metrics.

SummaryMonitoring subsurface fluid transport processes using electrical resistivity tomography (ERT) gained attention due to its sensitivity to variations in fluid content and temperature. Capturing these dynamic processes, such as tracer transport or pollutant dispersion, requires not only advanced imaging techniques but also efficient data acquisition strategies to ensure high-resolution and cost-effective monitoring. This study advances Optimal Experimental Design (OED) for transient flow monitoring using ERT. We developed three different strategies based on the Compare-R algorithm, and present as well as compare their results. By optimizing the selection of measurement configurations, these strategies aim to maximize the information content of the acquired data while minimizing redundancy and acquisition costs. The data-driven approach uses data from prior time steps to generate focusing masks. The model-driven approach incorporates transport models to target critical areas, and the hybrid approach combines both methods, iteratively refining the transport model based on the acquired data. We use synthetic studies to demonstrate the model-driven strategy”s advantage in identifying transport-affected areas, especially for rapid processes or longer monitoring intervals. The data-driven approach performs well when masked and process-affected areas overlap, but struggles with spatial delays otherwise. The hybrid OED effectively combines the strengths of both methods, detecting discrepancies between modeled and observed fluid flow, and refining transport simulations for subsequent steps, thereby ensuring reliable optimization for dynamic monitoring scenarios.

SummaryShear wave splitting, caused by seismic anisotropy, provides insights into convective processes in the mantle. While upper mantle anisotropy beneath the continents is well-resolved with high lateral resolution, its characterization beneath the oceans remains limited due to the paucity of seismic stations. In this study, we leverage the sensitivity of P waves that convert to S upon reflection at the surface (PS) to infer seismic anisotropy near their bounce point. We measure PS splitting automatically from a global dataset that includes all earthquakes with magnitudes ≥5.9 from January 1995 to present, collected from 25 datacenters, totaling approximately 5,800 events and 18 million three-component seismograms. After careful quality control, we obtain 889 PS splitting measurements from this dataset, mostly for regions without other splitting constraints. These regions include the South China Sea, the Philippine Sea, the western Pacific Ocean, the western Atlantic Ocean near the Caribbean, and the South Atlantic Ocean. Where independent SKS constraints are available, these are in general agreement with our PS results. Comparisons of PS fast directions to azimuthal seismic anisotropy derived from surface waves and the direction of absolute plate motion (APM) show a dominant contribution of asthenospheric deformation to the measurements. This agreement is particularly pronounced beneath the oceans, away from subduction zones, for instance, across much of the Pacific basin. In contrast, splitting patterns near subduction zones are generally more complex: fast directions are most often trench-parallel or oblique, and only rarely trench-perpendicular, suggesting that a trench-parallel component of shear deformation is commonly generated in subduction settings, as commonly inferred from SKS splitting. We make all measurements publicly available as a data product, along with detailed metadata to enable future work.

SummaryMass redistribution during earthquake rupture, along with subsequent wave propagation, perturbs Earth’s gravity field, generating so-called ‘prompt elasto-gravitational signals’ (PEGS) and ‘prompt gravity-strain signals’ (PGSS). These signals are detectable on accelerometers and gradiometers before P-wave arrivals, and therefore offer the potential for early-warning systems and rapid assessment of event magnitude and tsunamigenic risk. Despite their significance, numerical modelling of PEGS and PGSS has been restricted to 1D Earth models, assuming negligible effects of 3D heterogeneities.In this study, we utilise a Spectral-Infinite-Element method to compute these PEGS and PGSS for both 1D and 3D Earth models. Kernels of sensitivity to model heterogeneity are investigated using the adjoint method, which is re-formulated to write the self-gravitating adjoint equations directly in terms of density and elastic perturbations only. We find that PEGS and PGSS display the greatest sensitivity to model perturbations below the source and receiver; however, sensitivity is generally weak, validating previous assumptions that 1D models are sufficient.

SummaryThe strength of the subduction interface plays a major role in controlling subduction dynamics on both local and global scales. While previous studies have primarily examined interface strength in two-dimensional models, natural subduction zones are inherently three-dimensional, with interface strength varying along-strike due to spatial differences in factors such as sediment input. Here, we use the geodynamic code ASPECT to conduct fully dynamic 3D subduction models in which interface strength varies along-strike. We find that the interaction between strong and weak segments of the interface leads to a narrower range of convergence velocities while broadening the range of viable interface stresses compared to 2D or homogeneous 3D models. Stronger segments, when adjacent to weaker ones, exhibit increased convergence velocities. This promotes higher interface stresses and facilitates the subduction of otherwise stagnant strong segments. We find that the interface viscosity of the strong segment controls the baseline stress, whereas the viscosity contrast along-strike controls the magnitude of amplification of the stress due to velocity increases. The elevated interface stresses at strong segments also generate greater compressional forces in the overriding plate than expected from 2D models. Combined with along-strike variations in convergence velocity, this results in trench migration, with stronger segments displaying more advanced trench positions relative to weaker segments. Possible natural analogs include the Bolivian Orocline in the central Andes and the Lesser Antilles, both of which show enhanced overriding plate compression and trench advance in areas of reduced sediment supply.

SummaryRealistic models of earthquake sequences can be simulated by assuming faults governed by rate-and-state friction embedded in an elastic medium. Exploring the possibility of using such models for earthquake forecasting is challenging due to the difficulty of integrating Partial Differential Equation (PDE) models with sparse, low-resolution observational data. This paper presents a machine-learning-based reduced-order model (ROM) for earthquake sequences that addresses this limitation. The proposed ROM captures the slow/fast chaotic dynamics of earthquake sequences using a low-dimensional representation, enabling computational efficiency and robustness to high-frequency noise in observational data. The ROM’s efficiency facilitates effective data assimilation using the Ensemble Kalman Filter (EnKF), even with low-resolution, noisy observations. Results demonstrate the ROM’s ability to replicate key scaling properties of the sequence -namely the magnitude-frequency, moment-duration, and moment-area relationships- and to estimate the distributions of fault slip rate and state variable, enabling predictions of large events in time and space with uncertainty quantification. These findings underscore the ROM’s potential for forecasting and for addressing challenges in inverse problems for nonlinear geophysical systems.

SummaryThis study introduces a novel inversion approach to resolve the layered anisotropic structure of the Earth’s crust and upper mantle using harmonic patterns observed in Ps receiver functions (RFs). Designed for dense seismic networks, our method effectively captures the complexity and deformation of subsurface structures by leveraging the decomposition of RF patterns into five harmonic functions representing distinct terms in a harmonic regression model. Our approach combines residual weighting at individual and coherency weighting between stations to reduce noise and enhance structural signal, followed by a direct inversion of harmonic patterns alongside a multi-phase stack. This process is optimized using a Markov chain Monte Carlo framework to iteratively explore model space and define posterior distributions, which yields robust model parameter estimates with quantified uncertainties. We validated our approach with both synthetic models and real data from the dense SEISConn seismic transect in Northern Connecticut. The synthetic test highlighted the reliability of coherency weighting in reducing noise. Real data analysis revealed anisotropic and structural features consistent with geological expectations and previous studies, including shallow anisotropy in western Connecticut, a shallowing Moho towards the Hartford basin, and indications for a west-dipping layer beneath the Moho in the lithospheric mantle. Our approach offers a promising tool for revealing the details of anisotropic features beneath dense seismic profiles, facilitating insights into tectonic history and lithospheric deformation.

SummaryLarge-scale geodetic monitoring of volcanic earthquakes is essential for understanding the physical mechanisms governing volcanic activity and magma migration. Recently, a significant earthquake swarm occurred around the Scotia plate. Geodetic data of 14 permanent GNSS stations on the Antarctic Peninsula, King George Island, the South Sandwich Islands, and South America were collected and analyzed, to monitor the crustal deformation of King George Island, the expansion of Bransfield Strait, the drift of Scotia plate, and sea level anomalies. Tidal data from four permanent tide stations were analyzed to monitor sea level anomalies. Results showed that after earthquakes the King George Island’s movement speed increased tenfold and its direction altered by 90 degrees. Land surface fluctuations in southeast King George Island were observed a year before the earthquakes, followed by continuous uplift. A combinatorial model including a point pressure source and expanding dike fit well with new geodetic monitoring data, revealing the impact of volcanic activity on this region. Geodetic monitoring and modeling quantitatively depicted the pre-seismic, co-seismic, and post-seismic phases of geological changes, providing new evidence and insights into the complex geological structures.

SummaryA 220 km long active seismic profile crossing the Saronic and Corinth Gulfs was performed using 35 4C Ocean Bottom Seismographs (OBS) and 4 3C stand-alone land stations. We recorded shots fired along line at every 120 m from a 48-l airgun array of 51 bar-m power. The velocity model was developed by first break tomographic inversion, followed by kinematic and dynamic ray tracing. The final velocity model was used to prestack depth migrate the Common Receiver Gathers. Moho depth below the central Corinth Basin is located at 32 km, thinning to 22 km below the Lechaio Gulf, at the transition to the Saronic Basin. In the Saronic domain Moho is found at 19 to 21 km depth, whereas below we identified a low velocity upper mantle of Vp 7.5 to 7.6 km/s, extending to 34 km depth. This is a low velocity asthenosphere wedge that intruded from the Cyclades region below the Saronic Gulf, driving the volcanic activity. It does not extend in the Corinth Rift domain to the west, where Pn is 8.0 km/s. Two major extensional detachments were mapped along the profile: one to the NW in the central Corinth Gulf, east of Galaxidi, corresponding to the onshore Itea-Amfissa Detachment, having a throw of more than 2000 m; the other to the SE, west of Agios Georgios Island, separates the Cycladic Metamorphic Core Complex from the non-metamorphic internal nappes of the Hellenides. Thrusting to the NW observed in the SE Saronic Basin has doubled the thickness of the high velocity limestones. At 20 km SE of Agios Georgios Island a major dextral strike slip fault was mapped. Normal faults with more than 1 km throw are observed around the Isthmus of Corinth trending E-W and WNW-ESE. Maximum thickness of 2000 m of Middle Miocene-Quaternary sediments (Vp 2.1–2.8 km/s and 3.4 km/s) is observed in the Corinth Gulf and minimum of 200 m of Quaternary sediments above the western Cyclades. Beneath the volcanoes of Aegina, Methana, Paphsanias, and Sousaki the crust has a lateral Vp velocity increase of nearly 3 per cent, whereas depth migrated data show high reflectivity, indicating magma intrusions through the crust. Mesozoic limestones with Vp 5.2 to 5.8 km/s and 1800 to 2000 m thickness occur along the profile, corresponding probably to the Tripolis external carbonate platform, whereas limestones and flysch with Vp 4.3 to 5.5 km/s are overlying, probably corresponding to pelagic sequences like the Pindos nappe. The Saronic domain is characterized by arc-parallel NW-SE structures, hosting the volcanic arc, and is tectonically controlled by magma uplift and thermally triggered deformation. The Corinth domain is affected by E-W transverse-oblique structures of a rapidly developing continental rift, and deformation driven by fracturing of a brittle crust of the Central Hellenic Shear Zone.

SummaryTo better understand Eocene ophiolite emplacement and subduction initiation in northeastern Zealandia, we analysed ambient noise to image shallow (0–3 km) shear wave velocity structures of and beneath an ophiolite nappe in southern Grande Terre, New Caledonia. We assessed the uncertainties of each dispersion curve to obtain stable dispersion curves at short periods (<1 s) from a network of 17 seismic stations, whose average interstation distance is ∼15 km. We obtained 1D velocity profiles and interpolated them to generate 2D transects with a lateral resolution <2 km. Two velocity discontinuities were imaged at depths of 100 m and 400–700 m, representing surface regolith and the base of the ophiolite nappe, respectively. The ophiolite nappe is underlain by continental basement rocks in the centre of the island and sedimentary rocks near the east and west coasts. The base of the nappe shallows at ∼2.5° westward to the surface at its southwest flank. Based on the geometry of the ophiolite nappe, we suggest a down-going gravity-driven emplacement mechanism, and note similarities to allochthons in Reinga Basin and Raukumara Basin of northern New Zealand. The ophiolite nappe and underlying bedrock are more fractured on their east flank due to syn/post emplacement deformation, isostatic adjustment and present flexural bending.

SummarySeismological inversion traditionally targets either source parameters, such as location and moment tensor, or structural parameters, such as velocity and anisotropy. However, the natural formulation of Full-Waveform Inversion, often used for high-resolution structural model estimation, is to jointly invert for source and structural parameters. The common practice of holding source parameters, after initial estimation, fixed throughout the inversion inherently leads to biased solutions of the structural model, and vice versa. Whereas a joint inversion suffers from severe non-uniqueness, we demonstrate that leveraging the large amounts of data available from Distributed Acoustic Sensing (DAS) can yield robust and unbiased estimations of source and structural parameters, provided an appropriate misfit function and optimisation scheme are used. We show how the size of the data space and eventual convergence can be improved by supplementing the phase misfit objective function with amplitude information. To this end, we formulate a new misfit function, the normalised envelope. To support native DAS data implementations, we calculate the adjoint sources for the new misfit function when defined directly on strain or strain-rate data. We also show how a new approach to preconditioning as part of the L-BFGS optimisation scheme allows for effective updates of all parameters in the same iteration, despite enormous differences in their relative importance. We test our approach in a challenging synthetic noisy 2D scenario, showing a considerable reduction in source parameter errors and an improved S-wave velocity model. We also show a 3D synthetic case with an idealised DAS recording array, demonstrating a significant reduction of source parameter errors using realistic initial estimates and structural model errors. We argue that the proposed methodology can be used to improve the quality of earthquake catalogues and high-resolution structural models in seismically active regions, especially at the local-to-regional scale. None the less, computational cost remains a major challenge of the method.

SummaryTo understand the melt source of hotlines with asynchronous volcanoes, we investigate the lithospheric structure of the Cameroon Volcanic Line (CVL), an intraplate hotline without age progression stretching from the Atlantic Ocean into Central Africa. We analyze Bouguer gravity anomalies from the World Gravity Model 2012 using the 2‐D power spectrum techniques and 2-D forward modeling to estimate the crustal and lithospheric thickness. We find: (1) thin crust (20–30 km) beneath the oceanic CVL; (2) thick crust (30–43 km) beneath the continental CVL and the Oubanguides Belt, and thicker crust (43–50 km) beneath the Congo Craton; (3) thin lithosphere (90–120 km) beneath the oceanic CVL and thinner lithosphere (75–90 km) beneath the continental CVL; and (4) thicker lithosphere (150–234 km) beneath the Congo Craton. Our seismically constrained forward models reveal a delaminated body beneath the continental CVL and a sharp transition from thick lithosphere beneath the Congo Craton to thin lithosphere beneath the Oubanguides Belt. We interpret that the thin lithosphere beneath the continental CVL is a result of lithospheric delamination. The delaminated body in the uppermost mantle deflects rising mantle plume material, resulting in the Y-shaped distribution of continental volcanoes. Edge-Driven Convection (EDC) resulting from the sharp gradient in lithospheric thickness between the Congo Craton and the Oubanguides Belt focuses the plume material beneath thin lithosphere, producing the continental CVL. The southern volcanoes of the continental CVL are formed from the southward deflection of plume material by the delaminated body, with melt ascent facilitated by the lithospheric-scale Central African Shear Zone. The northward-directed plume material forms the distinct Biu Plateau, and the eastward-deflected plume material forms the Adamawa Plateau. With a continuous influx of plume material beneath the thin continental lithosphere, for mass to be conserved, part of the plume material defiles the gradient of the thicker oceanic lithosphere adjacent to the Congo Craton to flow oceanward. The oceanward flow of plume material is modulated by upwellings from EDC, producing the oceanic CVL, which explains the oceanward decrease in the timing of the onset of volcanism. We therefore conclude that only the continental CVL lacks age progression resulting from the complex interaction of the rising plume with the delaminated body and the lithospheric architecture.

SummaryInferring the spatio-temporal distribution of slip during earthquakes remains a significant challenge due to the high dimensionality and ill-posed nature of the inverse problem. As a result, finite-source inversions typically rely on simplified assumptions. Moreover, in the absence of ground-truth measurements, the performance of inversion methods can only be evaluated through synthetic tests. Laboratory earthquakes offer a valuable alternative by providing “simulated real data” and ground truth observations under controlled conditions, enabling a more reliable evaluation of source inversion procedures. In this study, we present static and quasi-static slip inversion results from data recorded during laboratory earthquakes. Each event is instrumented with 20 accelerometers along the fault, and the recorded acceleration data are used to invert for the slip history. We consider two different types of Green’s functions (GF): simplistic GF assuming a homogeneous elastic half-space and realistic GF computed by finite element modeling of the experimental setup. The inversion results are then compared to direct observations of fault slip and rupture velocity obtained independently during the experiments. Our results show that, regardless of the GF used, the inversions fit well with the data and result in small formal uncertainties of model parameters. However, only the inversion with realistic GF yields slip distributions consistent with the true fault slip measurements and successfully recovers the distribution of rupture velocity along the fault. These findings emphasize the critical role of GF selection in accurately resolving slip dynamics and highlight an important distinction in Bayesian inversion: while posterior uncertainty quantification is essential, it does not guarantee accuracy, especially if forward modeling uncertainties are not properly accounted for. Thus, confidence in inversion results must be paired with careful modeling choices to ensure physical reliability.

SummaryThe regional seismic travel-time (RSTT) model predicts travel times of regional seismic phases accounting for three-dimensional structure of the crust and the upper mantle on a global scale. Previous versions of the RSTT model have been implemented using nodes separated by ∼1° spacing across the globe. A regional-scale study using regional Pn and Pg travel times across Israel and the Middle East demonstrated that data driven, systematic grid refinement reduces travel time residuals and enhances resolution of smaller tectonic features in regions having dense ray coverage. High density Pn ray coverage in the western US, Europe, Middle East, and East Asia can likewise provide the resolution that allows systematic global grid refinement of the RSTT model. In this study, we use a large number of Pn ray paths originating from events located with an epicentral location uncertainty of 25 km (GT25) or better. We conduct targeted grid refinements at 1.0°, 0.5°, 0.25°, and 0.125° on a global scale, producing a refined RSTT model that yields a 21.6% reduction in median event location error in Europe and the Middle East, when compared with the original global RSTT model presented in Begnaud et al. (2021a). The new model also resolves finer tectonic structures in regions with high Pn ray density.

SummaryWe present and apply a pseudo trans-dimensional inversion method for 3D geometrical gravity inversion, in which the number of rock units, their geometry, and their density can vary during sampling. The method is designed for efficient exploration of the model space and to infer the presence and properties of units not directly observable but detectable with geophysical data. Sampling relies on a non-reversible Metropolis-Hastings algorithm, during which rock units can be added or removed from the model, interface geometries are perturbed using random fields, and densities are sampled from distributions informed by prior information. To visualise the space of sampled models and to aid interpretation, a workflow is proposed that combines dimensionality reduction with the clustering of models in families. The capabilities of the inversion method are evaluated using two synthetic cases. The first is a motivating example aimed at recovering an intrusion missing from the prior model. It features a horizontal layer-cake where fixed-dimensional inversion fails to adequately fit the data and sample models close to the true model, while the proposed pseudo trans-dimensional approach is much more successful. The second case investigates the recovery of two missing units and the capability to overcome prior model biases. Results show the potential of our method to infer the presence of unseen geological features such as intrusions. However, they suggest that with biased prior geological modelling, it may be challenging to infer with certainty the presence of more than two previously unknown rock units at depth.

SummaryWhen highly viscous fluids are present in a rock medium, the viscous effect of such fluids cannot be neglected in the propagation of elastic waves. In this paper, a fractional thermoporoelastic theory is newly proposed, which is a further improvement of the two-temperature generalized thermoporoelastic theory. Firstly, by introducing the Kelvin-Voigt model into the stress-strain constitutive equation, the viscous effect of highly viscous fluids is considered. Then, fractional derivatives are introduced into the heat conduction equations of the solid and fluid phase to consider the anomalous heat conduction caused by the viscous effect in rock media. Plane wave analysis method is adopted to obtain the phase velocity and attenuation factor of four longitudinal waves (P1, P2, T1, T2). Numerical results show that the introduction of fluid viscosity leads to the appearance of new relaxation peaks in the P wave at high frequencies, and the introduction of fractional derivatives causes a decrease in the phase velocity and attenuation factor of T waves. The results provide a reference for further research on the wave propagation in rock media containing highly viscous fluids.

SummaryInduced polarization (IP) effects in airborne electromagnetic (AEM) surveys have commonly been investigated in helicopter-borne systems, leaving both a bibliographic and application gap for fixed-wing configurations. This gap partly reflects the large relative number of helicopter compared to fixed wing AEM systems, but also the geometric complexity of fixed wing platforms. In these platforms, nine geometric parameters come into play: the pitch, roll, and yaw of both transmitter and receiver, plus the three-axis offsets between the coils. Shifts in these factors can distort the measured data in ways that aren’t uniquely attributable, making it hard to pinpoint whether negative recordings truly arise from IP or from geometry-related effects. The non-fixed geometry also complicates removal of the primary field, often requiring iterative processing steps that may suppress or alter spectral content linked to IP. With advances in airborne IP understanding from helicopter-borne systems, revisiting fixed-wing platforms is both timely and necessary. Part A of this two-part study addresses this issue using the TEMPEST™ fixed-wing system connecting numerical modelling with field evidence. A suite of synthetic two-layer models with variable resistivity and chargeability parameters was developed to evaluate the system’s sensitivity to polarizable structures. The experiments demonstrate that IP effects, including negative secondary field responses, can be reliably detected in fixed-wing AEM data, both in X and Z magnetic field components. The capacity of these systems to detect IP phenomena is, however, strongly dependent on the electrical conductance of the environment. For instance, both fixed-wing and helicopter-borne systems, elevated near-surface conductance enhances the amplitude of purely electromagnetic induction currents, which in turn can dominate the recorded response and obscure the comparatively weaker polarization currents. More in general, IP detectability depends on the strength of the EM response generated by induction currents flowing elsewhere, which can dominate the small reverse current flow from a polarizable target. This highlights the critical role of near-surface conductivity in controlling the expression of IP responses and underscores the need to carefully account for these factors when interpreting survey data. The synthetic results are then connected with field-scale observations from a subset of the AusAEM dataset, over 470 000 line-km of TEMPESTTM data, where negative responses align with areas of low shallow conductance, confirming the simulation results. These finding open the way to the Part B of this study, where TEMPESTTM data are inverted taking into account IP and compared with helicopter-borne results and geological information.

SummaryAcoustic signals can couple to the ground, providing an opportunity to use seismic stations to investigate airborne sources. The study of Bishop et al. (2022) used wavefield simulations in a fluid-solid medium to quantify the role of topography on the seismic (ground) recordings of a monopole source in the air. We build upon this study by linking wavefield forward modeling with the source estimation code MTUQ, which can accommodate point forces or moment tensors in a solid medium, as well as sources in the air (or water) if they are enabled by the forward-modeling solver. We perform a series of synthetic numerical experiments to demonstrate that a dipole airborne source can be estimated using ground-based receivers, even within the presence of realistic topography. We investigate the influence of receiver coverage, topography, and assumed source location on the estimated results. The established capabilities raise the prospects for future efforts to estimate dipole sources in 3D models that include heterogeneity in the air and the earth in addition to topography.

AbstractSeismic simulation is fundamental for understanding earthquake physics and mitigating seismic hazards, but accurate seismic modeling requires fine computational grids, imposing severe memory and computational challenges. Traditional modeling solvers, relying on single-precision floating-point 32-bit (FP32) and scalar register-based computation, suffer from excessive memory consumption, low memory access efficiency, and limited computational efficiency. Compared with FP32, half-precision floating-point 16-bit (FP16) reduces memory consumption by 50% and improves memory access efficiency; relative to scalar registers, ARM’s Scalable Vector Extension (SVE) registers provide vectorized single-instruction multiple-data (SIMD) capabilities, significantly accelerating computation. However, leveraging the advantages of FP16 and SVE involves challenges such as FP16 overflow/underflow, SVE stencil adaptation, and SVE data misalignment from FP16 storage with FP32 arithmetic. Therefore, this study proposes three approaches on the ARM architecture: FP16-based, SVE-accelerated, and FP16-SVE hybrid; each is designed to tackle the respective challenges while exploiting FP16 memory efficiency and SVE computational acceleration. Correspondingly, the three solvers are implemented, validated, and benchmarked on both hypothetical models and real-world earthquake scenarios. The results of these solvers show near-perfect agreement with the reference solver, confirming their accuracy across diverse seismic scenarios. Moreover, the FP16-SVE hybrid solver halved memory usage and achieved up to 3× computational speedup, delivering more than 2.3× acceleration in the real-world earthquake simulation. The gains in high efficiency of memory and computation highlight the capability of the FP16-SVE hybrid solver to support large-scale, real-time seismic simulations and efficient earthquake hazard assessment on ARM platforms.

AbstractSurface waves such as Rayleigh, Love and Scholte waves can exhibit dispersion, i.e., variations in phase velocity with wavelength as a function of frequency. This property enables the inversion of 1D models of seismic velocity and density in the subsurface. Conventional deterministic and stochastic inversion schemes are widely applied to surface wave data but face two main challenges. The first is the identification of dispersion curves for fundamental and higher modes on wavefield-transformed images, which is often done manually. The second is the quantification of uncertainty, which can be computationally expensive in stochastic approaches or limited to data-propagated uncertainty in deterministic inversions. Our objectives are to (1) eliminate the need for manual or automatic dispersion curve picking, and (2) directly infer ensembles of 1D velocity models - and their associated uncertainties - from the full velocity spectrum, i.e., the complete dispersion image containing all modes. To this end, we employ Bayesian Evidential Learning, a predictive framework that reproduces experimental data from prior information while allowing prior falsification. In our application, ensembles of prior Earth models are sampled to predict 1D subsurface structures in terms of seismic velocity and, where applicable, attenuation from near-surface seismic wave data. This approach bypasses traditional inversion schemes and provides a computationally efficient tool for uncertainty quantification.