Seismic Imaging Research Articles

Estimation of Water Velocity and Receiver Position Changes in Time-Lapse Seismic Using Machine Learning

Non-repeatability is one of the main obstacles in time-lapse seismic, as it significantly degrades the interpretation of reservoir-related signals. Correcting the data variations caused by non-repeatability (4D noise) is of paramount importance, which usually requires the estimation of the changing parameter. ln this paper, we propose a Machine Learning (ML) workflow for the quantitative estimation of two types of 4D noise: changes in the speed of sound in water and receiver lateral positions. A synthetic database, modeled from a velocity model estimated from the Brazilian pre-salt, and containing many time-lapse seismic surveys was generated for the supervised training of ML models. lnput samples consist of subsets of common-shot seismograms. We studied many combinations of ML regression algorithms and feature extraction techniques, for scenarios where data is contaminated, or not, with Gaussian random noise. The regression algorithms considered were Fully-Connected Neural Network, Extreme Gradient Boosting, and Bayesian Ridge. Four rectangular crops of the common-shot seismogram were tested as input features: full seismogram, first half of time samples, 11 smallest-offset traces, and the time samples focusing on the first arrivals region. The combination with the best trade-off between accuracy and model complexity is the Bayesian Ridge fed with the 11 smallest-offset traces, which estimated position and velocity time-lapse changes with median accuracy of 0.115 m and 0.017 m/s for the case with Gaussian noise. Besides the correction of repeatability-related variations, our results are useful in the 4D Full-Waveform lnversion which needs accurate parameters to produce good seismic images.

Brief communication: On the potential of seismic polarity reversal to identify a thin low-velocity layer above a high-velocity layer in ice-rich rock glaciers

Abstract. Seismic refraction tomography is a commonly used technique to characterise rock glaciers, as the boundary between unfrozen and ice-bearing layers represents a strong impedance contrast. In several rock glaciers, we observed a reversed polarity of the waves refracted by an extended ice-bearing layer compared to direct-wave arrivals. This phase change may be related to the presence of a thin low-velocity layer (LVL), such as fine- to coarse-grained sediments, above a thicker ice-rich layer. Our results are confirmed by the modelling and analysis of synthetic seismograms to demonstrate that the presence of a low-velocity layer can produce a polarity reversal on the seismic gather.

Waveform modelling of hydroacoustic teleseismic earthquake records from autonomous Mermaid floats

SUMMARY We present a computational technique to model hydroacoustic waveforms from teleseismic earthquakes recorded by mid-column Mermaid floats deployed in the Pacific, taking into consideration bathymetric effects that modify seismo-acoustic conversions at the ocean bottom and acoustic wave propagation in the ocean layer, including reverberations. Our approach couples axisymmetric spectral-element simulations performed for moment-tensor earthquakes in a 1-D solid Earth to a 2-D Cartesian fluid–solid coupled spectral-element simulation that captures the conversion from displacement to acoustic pressure at an ocean-bottom interface with accurate bathymetry. We applied our workflow to 1129 seismograms for 682 earthquakes from 16 Mermaids (short for Mobile Earthquake Recording in Marine Areas by Independent Divers) owned by Princeton University that were deployed in the Southern Pacific as part of the South Pacific Plume Imaging and Modeling (SPPIM) project. We compare the modelled synthetic waveforms to the observed records in individually selected frequency bands aimed at reducing local noise levels while maximizing earthquake-generated signal content. The modelled waveforms match the observations very well, with a median correlation coefficient of 0.72, and some as high as 0.95. We compare our correlation-based traveltime measurements to measurements made on the same data set determined by automated arrival-time picking and ray- traced traveltime predictions, with the aim of opening up the use of Mermaid records for global seismic tomography via full-waveform inversion.

Highlights on mantle deformation beneath the Western Alps with seismic anisotropy using CIFALPS2 data

Abstract. There are still open questions about the deep structure beneath the Western Alps. Seismic velocity tomographies show the European slab subducting beneath the Adria plate, but all these images did not clarify completely the possible presence of tears, slab windows, or detachments. Seismic anisotropy, considered an indicator of mantle deformation and studied using data recorded by dense networks, allows a better understanding of mantle flows in terms of location and orientation at depth. Using the large amount of shear wave-splitting and splitting-intensity measurements available in the Western Alps, collected through the CIFALPS2 temporary seismic network, together with already available data, some new patterns can be highlighted, and gaps left by previous studies can be filled. Instead of the typical seismic anisotropy pattern parallel to the entire arc of the Western Alps, this study supports the presence of a differential contribution along the belt that is only partly related to the European slab steepening. A nearly north–south anisotropy pattern beneath the external Western Alps, a direction that cuts the morphological features of the belt, is clearly found with the new CIFALPS2 measurements. It is, however, confirmed that the asthenospheric flow from central France towards the Tyrrhenian Sea is turning around the southern tip of the European slab.

Reconstruction and denoising of high-dimensional seismic data via Frobenius-nuclear mixed norm constraints

Abstract The seismic data acquisition design with ‘two-wide and one-high’ geometry effectively improves the imaging quality of seismic records. However, when data are acquired in the field, complex near surface conditions and environmental factors can introduce a variety of noises and gaps in seismic data, impacting the accuracy of seismic imaging. Currently, the method of low-rank matrix/tensor completion is commonly employed for data reconstruction after normal moveout (NMO). In a complex subsurface medium, common midpoint data processed with NMO may not satisfy the linear or quasi-linear assumptions within local data windows. Therefore, this paper exploits the inherent low-rank structure of high-dimensional data to propose a high-dimensional tensor completion method under the Frobenius-nuclear mixed norm constraint (FN-TC). This method unfolds the 4D data tensor into the frequency-space domain along its modes (m, n) and subsequently imposes a non-convex Frobenius-nuclear mixed norm constraint on the unfolded approximate matrices. This approach closely approximates the rank function of the factor matrices, thereby enhancing the accuracy of data modeling. Theoretical and practical studies demonstrate that the novel FN-TC approach can effectively reconstruct high-dimensional seismic data and suppress noise, thereby providing data support for subsequent high-precision seismic imaging.

Refining tomography with generative neural networks trained from geodynamics

SUMMARY Inverse problems occur in many fields of geophysics, wherein surface observations are used to infer the internal structure of the Earth. Given the non-linearity and non-uniqueness inherent in these problems, a standard strategy is to incorporate a priori information regarding the unknown model. Sometimes a solution is obtained by imposing that the inverted model remains close to a reference model and with smooth lateral variations (e.g. a correlation length or a minimal wavelength are imposed). This approach forbids the presence of strong gradients or discontinuities in the recovered model. Admittedly, discontinuities, such as interfaces between layers, or shapes of geological provinces or of geological objects such as slabs can be a priori imposed or even suggested by the data themselves. This is however limited to a small set of possible constraints. For example, it would be very challenging and computationally expensive to perform a tomographic inversion where the subducting slabs would have possible top discontinuities with unknown shapes. The problem seems formidable because one cannot even imagine how to sample the prior space: is each specific slab continuous or broken into different portions having their own interfaces? No continuous set of parameters seems to describe all the possible interfaces that we could consider. To circumvent these questions, we propose to train a Generative Adversarial neural Network (GAN) to generate models from a geologically plausible prior distribution obtained from geodynamic simulations. In a Bayesian framework, a Markov chain Monte Carlo algorithm is used to sample the low-dimensional model space depicting the ensemble of potential geological models. This enables the integration of intricate a priori information, parametrized within a low-dimensional model space conducive to efficient sampling. The application of this approach is demonstrated in the context of a downscaling problem, where the objective is to infer small-scale geological structures from a smooth seismic tomographic image.

Integrated geological and geophysical workflow for structural modelling: a case study from the contractional foothills zone of the Colombian Eastern Cordillera

In fold-and-thrust belts where there is a high degree of structural complexity, artificial geometrical distortions are often present on seismic reflection profiles. These need to be minimized during modelling. We document a workflow in which depth mapping, velocity model building, well calibration and cross-section balancing are integrated into the seismic interpretation process to generate trustworthy structural models in complex zones. The proposed methodology is exemplified by a case study from the foothills zone of the Colombian Eastern Cordillera. In addition, sequential kinematic restoration of the modelled structure allowed an evaluation to be made of hydrocarbon migration routes during the period between the Oligocene and the middle Miocene. Following the previously mentioned workflow, we document a failed exploratory case study where all elements of the petroleum system are present except the trap. In this context, the documented case is a typical velocity pull-up. From this and published case studies, we conclude that in Andean settings and probably most on land contractional-foothill settings, the use of seismic images only does not provide enough evidence for the presence of traps and additional surface geological signatures must be documented. The proposed workflow therefore appears to be a useful tool for evaluating the exploration risk in structurally complex fold-and-thrust belt settings.

Seismic imaging of the complex geological structures in the southwestern edge of the Western limb, Bushveld Complex through focusing pre‐stack depth migration of legacy 2D seismic data

Seismic imaging of the complex geological structures in the southwestern edge of the Western limb, Bushveld Complex through focusing pre‐stack depth migration of legacy 2D seismic data

Characteristics of the Fault Damage Zone From High‐Resolution Seismic Imaging Along the Palos Verdes Fault, California

AbstractThe distribution and intensity of fault damage zones provides insight into fault activity and its relationship to fluid flow in the crust. Presently, measures of the in‐situ distribution of fault damage remain limited and along‐strike studies are rare. This study focuses on an offshore section Palos Verdes Fault damage zone that spans 28 km, near Los Angeles, California. To investigate the previously unresolved shallow (∼400 m below the seafloor) fault damage zone we use densely spaced (∼500 m line separation) newly collected sparker multichannel seismic lines and sub‐bottom profiles. The combination of high‐resolution acquisition methods and specialized seismic processing workflows provide improved imaging of shallow faulting. We apply a multi‐trace similarity technique to identify discontinuities in the seismic data that may be attributed to faults and fractures. This fault detection approach reveals diverse fault damage patterns on adjacent seismic profiles. However, a discernible damage zone pattern emerges by stacking multiple damage detection profiles along strike. We find that peak damage identified in this way corresponds to the active main fault strand, confirmed in this study, and thus the technique may be useful for identifying active fault strands elsewhere. Additionally, we observe that the variable width of the damage zone along strike is controlled by fault obliquity. Furthermore, our observations reveal a correlation between fault damage and seafloor fluid seeps visible in the water column, suggesting that damage plays a role in controlling fluid flow around the fault.

Constraining Displacement Magnitude on Crustal‐Scale Extensional Faults Using Thermochronology Combined With Flexural‐Kinematic and Thermal‐Kinematic Modeling: An Example From the Teton Fault, Wyoming, USA

AbstractConstraining the geometry and displacement of crustal‐scale normal faults has historically been challenging, owing to difficulties with geophysical imaging and inability to identify precise cut‐offs at depth. Using a modified workflow previously applied to contractional systems, flexural‐kinematic (Move) and thermal‐kinematic (Pecube) models are integrated with apatite (U‐Th)/He (AHe) and apatite fission track (AFT) data from Teton footwall transects to constrain total Teton fault displacement (Dmax). Models with slip onset at ∼10 Ma and flexure parameters that best match the observed Teton flexural profile require Dmax > 8 km to produce young (<10 Ma) AHe ages observed at low elevation footwall positions in the Tetons. For the same slip onset, models with Dmax of 11–13 km provide the best match to observed AHe data, but displacements ≥16 km are required to produce observed AFT ages (13.6–12.0 Ma) at low elevations. A more complex model with slow slip onset at ∼25 Ma followed by faster slip at ∼10 Ma yields a good match between modeled and observed AHe ages at a Dmax of 13–15 km. However, this model predicts low elevation AFT ages 6–8 Ma older than observed ages, even at Dmax values of 16–17 km. Based on this analysis and integration with previous studies, we propose a unified evolution wherein the Teton fault likely experienced 11–13 km of Miocene‐recent displacement, with AFT data likely indicating a pre‐to early Miocene cooling history. Importantly, this study highlights the utility of using integrated flexural‐ and thermal‐kinematic models to resolve displacement histories in extensional systems.

Geophysical and numerical approaches to solving the mechanisms of landslides triggered by earthquakes: A case study of Kahramanmaraş (6 February, 2023)

Many landslides occurred in the 11 Turkish provinces affected by the earthquakes centered in Kahramanmaraş in February 2023, and many people lost their lives due to those landslides. To reduce the risks and hazards brought on by landslides induced by earthquakes, it is essential to identify the mechanisms causing landslides to arise in such situations. In this study, the geological strength index and both seismic refraction tomography and electrical resistivity tomography (ERT) geophysical methods were used to determine the failure mechanism of layered sedimentary units defined as flysch (claystone-siltstone-sandstone-marl). The shear strength reduction finite element method (SSR-FEM) was used with the RS2 computer program for model solutions. One of the geophysical methods, ERT, facilitated an interpretation of the cyclic failure mechanism that occurs in flysch, and it was found to be congruent with numerical models. Landslide slip circle depth, water-saturated layers, and resistivity values were determined with the ERT method. Seismic methods were employed to determine the velocity values (Vp, Vs) of the landslip mass. These velocity data, along with the elasticity modulus and Poisson ratios derived from them, were utilized in our numerical model. In the maximum shear strength and horizontal displacement analyses performed on the section line determined before the landslide, the natural state and the situation under the influence of groundwater were examined and the strength reduction factor (SRF) was determined as 1.35 and 1.24, respectively. In an evaluation conducted under the influence of the maximum horizontal ground acceleration (0.22 g) determined by taking into account the seismicity of the study area, it was concluded that instability occurred in the flysches and the SRF value was 0.72. An important finding of this study is that the results of geophysical methods used to determine the lithological properties of sedimentary rock masses and understand the causes of landslides align closely with those of numerical analysis methods, providing high accuracy. These methods also support new approaches, particularly in identifying landslip boundaries and detecting areas with high water saturation. It has been discovered that earthquakes actively contribute to the instability of these kinds of soils.

Multi‐Scale Seismic Imaging of the Ridgecrest, CA, Region With Waveform Inversion of Regional and Dense Array Data

AbstractWe develop an inversion procedure for deriving multi‐scale velocity models with waveform inversions of earthquake and ambient noise data at multi‐frequency bands recorded by regional and dense sensor configurations. The method is applied for the area around the 2019 Ridgecrest earthquake rupture zones, utilizing data recorded by regional stations and dense 2D and 1D arrays with station spacings of ∼5 km and ∼100 m, respectively. Starting with regional Vp, Vs models and locations of Ridgecrest aftershocks, the velocity models and event locations are improved iteratively by inversions of waveforms recorded by regional stations and the 2D array, using a minimum spectral element size of ∼600 m. Waveforms from local events recorded by dense 1D arrays across the M7.1 rupture zone with frequencies of up to 10 Hz are used to resolve small‐scale features of the rupture zone and shallow crust with a local spectral element size of 80 m. The refined models provide self‐consistent descriptions of the rupture zone and the shallow crust embedded in the regional structures. The results reveal pronounced low Vs and high Vp/Vs in the M6.4 and M7.1 rupture zones coinciding with concentrations of seismicity, and also around the Garlock fault and in several local basins. We also observe clear velocity contrasts across the Garlock fault with polarity reversals along strike and with depth. The obtained multi‐scale velocity models can be used to improve derivations of earthquake source properties, simulations of dynamic ruptures and ground motions, and the understanding of fault and tectonic processes in the region.

The Tectonic Structure and Evolution of the Potiguar‐Ceará Rifted Margin of Brazil

AbstractThe Brazilian Equatorial Margin (BEM) is interpreted as a transform margin, where the last segment opened during Gondwana rifting. However, margin evolution, and break‐up age remain unconstrained. We interpret >10k km of crustal‐scale seismic images extending along ∼600 km of the margin calibrated with drillholes. We determine the style and timing of tectonics across the rift system. We link changes in crustal‐scale structure and age of sediment deposits to interpret variations with the style of extension and intensity of thinning across the BEM. Observations support a rift evolution where deformation is initially distributed forming a shallow basin, subsequently focusses, and later migrates basin‐ward forming the deep‐water domain. We interpret that tectonic activity started ∼140–136 Ma and stopped earlier in the shallow basin causing minor thinning, than in the deep‐water domain with a ∼60 km wide area with 4–8 km thick crust extended in Late Aptian to Early Albian (116–110 Ma). Constraints from seismic and drilling help define an abrupt continent to ocean transition (COT) where continental crust may be abutted by oceanic crust, and breakup occurred at early Albian time. Basin sedimentation from the onset to the Late Aptian is continental, indicating an isolated environment disconnected from Atlantic oceans. During late‐most Aptian to Early Albian basin sedimentation changes and indicates a comparatively rapid marine water infill. Rifting of the BEM is not dominated by transcurrent deformation as previously inferred, with strike‐slip faulting limited to comparatively small sectors, whereas most of the margin extended by normal faulting deformation.

SeisResoDiff: Seismic resolution enhancement based on a diffusion model

High resolution of post-stack seismic data assists in better interpretation of subsurface structures as well as high accuracy of impedance inversion. Therefore, geophysicists consistently strive to acquire higher resolution seismic images in petroleum exploration. Although there have been successful applications of conventional signal processing and machine learning for post-stack seismic resolution enhancement, there is limited reference to the seismic applications of the recent emergence and rapid development of generative artificial intelligence. Hence, we propose to apply diffusion models, among the most popular generative models, to enhance seismic resolution. Specifically, we apply the classic diffusion model—denoising diffusion probabilistic model (DDPM), conditioned on the seismic data in low resolution, to reconstruct corresponding high-resolution images. Herein the entire scheme is referred to as SeisResoDiff. To provide a comprehensive and clear understanding of SeisResoDiff, we introduce the basic theories of diffusion models and detail the optimization objective's derivation with the aid of diagrams and algorithms. For implementation, we first propose a practical workflow to acquire abundant training data based on the generated pseudo-wells. Subsequently, we apply the trained model to both synthetic and field datasets, evaluating the results in three aspects: the appearance of seismic sections and slices in the time domain, frequency spectra, and comparisons with the synthetic data using real well-logging data at the well locations. The results demonstrate not only effective seismic resolution enhancement, but also additional denoising by the diffusion model. Experimental comparisons indicate that training the model on noisy data, which are more realistic, outperforms training on clean data. The proposed scheme demonstrates superiority over some conventional methods in high-resolution reconstruction and denoising ability, yielding more competitive results compared to our previous research.

High-precision seismic imaging for complex deep structures in the hydrocarbon exploration using a coherent-stacking-based least-squares migration

High-precision seismic imaging for complex deep structures in the hydrocarbon exploration using a coherent-stacking-based least-squares migration

Computation of Green’s Function in a Strongly Heterogeneous Medium Using the Lippmann–Schwinger Equation: A Generalized Successive Over-Relaxion plus Preconditioning Scheme

The computation of Green’s function is a basic and time-consuming task in realizing seismic imaging using integral operators because the function is the kernel of the integral operators and because every image point functions as the source point of Green’s function. If the perturbation theory is used, the problem of the computation of Green’s function can be transformed into one of solving the Lippmann–Schwinger (L–S) equation. However, if the velocity model under consideration has large scale and strong heterogeneity, solving the L–S equation may become difficult because only numerical or successive approximate (iterative) methods can be used in this case. In the literature, one of these methods is the generalized successive over-relaxation (GSOR) iterative method, which can effectively solve the L–S equation and obtain the desired convergent iterative series. However, the GSOR iterative method may encounter slow convergence when calculating the high-frequency Green’s function. In this paper, we propose a new scheme that utilizes the GSOR iterative with a precondition to solve the complex wavenumber L–S equation in a slightly attenuated medium. The complex wavenumber with imaginary components localizes the energy of the background Green’s function and reduces its singularity by enabling exponential decay. Introduction of the preconditioning operator can further improve the convergence speed of the GSOR iterative series. Then, we provide a preconditioned generalized successive over-relaxation (pre-GSOR) iterative format. Our numerical results show that if an appropriate damping factor and a proper preconditioning operator are selected, the method presented here outperforms the GSOR iterative for the real wavenumber L–S equation in terms of the convergence speed, accuracy, and adaptation to high frequencies.

Integrated Geophysical and Geotechnical Methods for Pre-Foundation Investigations at Lagos State University of Science and Technology, Ikorodu, Southwestern Nigeria

The recurrence of building collapses and resulting casualties is overwhelming in Lagos State as the state records the highest cases in recent times. This work used geophysical (Ground magnetic (GM) and Electrical resistivity imaging (ERI)) and geotechnical (Dynamic cone penetration test (DCPT)) methods in delineating the subsurface conditions for its competence in foundation structures. Ten profiles of 200 m ground magnetic data were collected at 2 m spacing likewise for Wenner arrays on each profile. DCPT data were collected at points 0 m, 50 m, 100 m,150 m, and 200m to the depth of 3 m on each of the profiles occupied by GM and ERI. The data obtained were processed and interpreted using RES2DINV software for ERI, Excel package for both Magnetic and DCPT. The ground magnetic results gives the minimum relative magnetic intensity as -1689.78 nT representing possible fault/fractured zone and maximum value of 2943 nTindicating area of competence for foundation structure. ERI results gives a resistivity value ranging from 2 Ωm – 7602 Ωm depicting four zones of incompetent/sub-competent clay/sandy clay, moderately competent clayey sand, competent sand and highly competent lateritic sand. DCPT results is in the average range 11.20% to 39.87% representing Clay materials fairly competent for foundation structures, sand and well graded gravel that is competent for construction purposes. To prevent building collapsed in the study area, areas of incompetent materials should be excavated, refill and compact with lateritic sand before foundation is put in place.

Present day mantle structure from global mantle convection models since the Cretaceous

SUMMARY Using forward mantle convection models starting at 140 Ma, and assimilating plate reconstructions as surface velocity boundary condition, we predict present-day mantle structure and compare them with tomography models, using geoid as an additional constraint. We explore a wide model parameter space, such as different values of Clapeyron slope and density change across 660 km, density and viscosity of the thermochemical piles at the core–mantle boundary (CMB), internal heat generation rate, and model initiation age. We also investigate the effects of different strengths of a weak layer below 660 km and weaker asthenosphere and slabs. Our results suggest that slab structures at different subduction zones are sensitive to the viscosity of the asthenosphere, strength of slabs, values of Clapeyron slope and the density and viscosity of the thermochemical piles, while different internal heat generation rates do not affect the slab structures. We find that with a moderately weak asthenosphere ($10^{20}$ Pa·s) and strong slabs, the predicted slab structures are consistent with the tomography models, and the observed geoid is also matched well. Moreover, our models successfully reproduce the degree-2 structure of the lower mantle beneath Africa and the Pacific, also known as Large Low Shear Velocity provinces (LLSVPs). A moderate Clapeyron slope of −2.5 MPa K−1 at 660 km aids in slab stagnation while higher values result in massive slab accumulation at that depth, ultimately leading to slab avalanches. We also find that the convective patterns in the thermal and thermochemical cases with slightly denser LLSVPs are similar, although the geoid amplitudes are lower for the latter. However, with more dense LLSVPs, the slabs cannot perturb them and no plumes are generated. Plumes arise as thermal instabilities from the edges of the LLSVPs, when cold and viscous slabs perturb them. While our predicted plume locations are consistent with the observed hotspot locations, matching the plume structures in tomography models is difficult. These plumes are essential in fitting the finer features of the observed geoid. In longer-duration models, more voluminous subducted material reaches the CMB, which tends to erode the LLSVPs significantly, and yields a poor fit to the observed geoid. Our results suggest that with the presence of a thin, moderately weak layer below 660 km, a slightly dense LLSVP, and Clapeyron slope of −2.5 MPa K−1, the velocity anomalies in seismic tomography and the long-wavelength geoid can be matched well. One of the limitations of our models is that the assimilated plate motion history may be too short to overcome arbitrary initial conditions effects. Also, assimilated true plate velocities in our models may not represent the true convective vigour of the Earth.

Efficient modeling of fractional Laplacian viscoacoustic wave equation with fractional finite-difference method

The fractional viscoacoustic/viscoelastic wave equation, which accurately quantifies the frequency-independent anelastic effects, has been the focus of seismic industry in recent years. The pseudo-spectral (PS) method stands as one of the most widely used numerical methods for solving the fractional wave equation. However, the PS method often suffers from low accuracy and efficiency, particularly when modeling wave propagation in heterogeneous media. To address these issues, we propose a novel and efficient fractional finite-difference (FD) method for solving the wave equation with fractional Laplacian operators. This method develops an arbitrary high-order FD operator via the generating function of our fractional FD (F-FD) scheme, enhancing accuracy with L2-optimal FD coefficients. Similar to classic FD methods, our F-FD method is characterized by straightforward programming and excellent 3D extensibility. It surpasses the PS method by eliminating the need for Fast Fourier Transform (FFT) and inverse-FFT (IFFT) operations at each time step, offering significant benefits for 3D applications. Consequently, the F-FD method proves more adept for wave-equation-based seismic data processes like imaging and inversion. Compared with existing F-FD methods, our approach uniquely approximates the entire fractional Laplacian operator and stands as a local numerical algorithm, with an adjustable F-FD operator order based on model parameters for enhanced practicality. Accuracy analyses confirm that our method matches the precision of the PS method with a correctly ordered F-FD operator. Numerical examples show that the proposed method has good applicability for complex models. Finally, we have carried out reverse time migration on the Marmousi-2 model, and the imaging profiles indicate that the proposed method can be effectively applied to seismic imaging, demonstrating good practicability.

The Fagradalsfjall and Sundhnúkur Fires of 2021–2024: A single magma reservoir under the Reykjanes Peninsula, Iceland?

AbstractThe Reykjanes Peninsula (RP) hosts several volcanic lineaments that have been periodically active over the last 4000 years. Since 2021, following a ca. 800‐year quiescence, eight eruptions have occurred on the RP, with more expected in the future. To better understand the origins of this renewed volcanism and help forecast future eruptions, we examine (i) if the ongoing volcanism is fed from a single or multiple magma storage zone(s) or from several smaller reservoirs and; (ii) where the zone(s) are located (i.e. mantle or lower or upper crustal depths). Using major and trace element geochemistry, oxygen isotopes, and seismic tomography we rule out a single, RP‐scale, deep‐seated magma storage zone. Instead we propose the presence of a ca. 10‐km‐wide region of crustal‐level (9–12 km) magma accumulation beneath the Fagradalsfjall volcanic lineament that fed both the 2021–23 eruptions of the Fagradalsfjall Fires and the 2023–24 eruptions of the Sundhnúkur Fires.