The 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, multiscattering 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.

Conventional reflection waveform inversion solves a two-parameter seismic inverse problem alternately for subsurface reflectivity and acoustic background velocity as the model parameters. It seeks to reconstruct a low-wavenumber velocity model of the subsurface from pure reflection data cyclically, through alternating migration and tomography loops, such that the remodelled data fits the observed data. Low-resolution seismic images with unpreserved amplitudes, full-wave inconsistency in the short-offset data and cycle skipping in the long-offset are perceived as the main reasons for suboptimal tomographic updates and slow convergence in conventional reflection waveform inversion. In the context of one-way reflection waveform inversion, this paper addresses the listed limitations through four main components. First, it augments one-way reflection waveform inversion with a computationally affordable preconditioned least-squares wave equation migration algorithm to ensure high-resolution reflectors with preserved amplitudes. Second, the paper verifies how well the full-wave consistency condition in the short-offset data is satisfied in one-way reflection waveform inversion and suggests muting inconsistent short-offset residual waveforms in the tomography loop to attenuate their adverse imprint. Third, the paper suggests extending the migration offset beyond short offsets to improve both the illumination and the signal-to-noise ratio of the reflectors. Fourth, the paper presents a data-selection algorithm to exclude the damaging effect of the cycle-skipped long-offset data in the tomography loop. The effectiveness of the proposed one-way reflection waveform inversion algorithm is finally validated through three numerical examples, demonstrating its capability to recover high-fidelity tomograms.

As seismic migration is increasingly applied to more and more complex media, more sophisticated imaging techniques are required to generate accurate images of the subsurface. Currently, the best results for imaging are achieved by least-squares migration methods, such as least-squares reverse time migration and full-wavefield migration (FWM). These methods iteratively update the image to minimize the misfit between the forward modelled wavefield and the recorded data at the surface. However, a key challenge for these techniques is the speed of convergence. To accelerate the speed of convergence, pre-conditioning is commonly applied. The most common pre-conditioner is the reciprocal of the Hessian operator. However, this operator is computationally expensive to calculate, making it difficult to apply directly. In this paper, we present a novel, alternative, pre-conditioner for FWM. This pre-conditioner is based on applying Galerkin projections to a linear system, which projects the system onto a set of known basis vectors. To find an appropriate set of basis vectors for this approach we apply proper orthogonal decomposition (POD) to a set of partial solutions of the linear system. The resulting method gives an approximation to the pseudo-inverse based on these basis vectors. To test this technique, which we name model-order reduced FWM (MOR-FWM), we apply it to the synthetic Marmousi model as well as to field data from the Vøring basin in Norway. For these examples, we show that MOR-FWM yields an improved data-misfit compared to the standard FWM approach. In addition, we show that the result for the field data case can be improved by normalizing the partial solutions before applying POD.

Imaging and inversion of land seismic data affected by complex weathering layers near the surface are challenging. When the data are additionally subsampled for economical reasons such as monitoring of sequestrated carbon dioxide and hydrogen, the problem is further exacerbated due to the combined influence of subsampling and weathering layers. First, interpolation performs poorly because the weathering layers reduce the data’s coherency. Second, near-surface corrections require knowledge of the subsurface model, separation between primaries and multiples, as well as subsurface velocity estimation, which are difficult to perform from subsampled data. To overcome these hurdles, we combine seismic interpolation and statics estimation into a joint single rank-reduction-based algorithm. To our knowledge, this is the first time that this has been done. Our method simultaneously accounts for the weathering and subsampling effects, which both contribute to the low-rank (LR) structure destruction typically associated with statics-free densely sampled data, to provide accurate reconstruction. Because an LR approximation is used for statics estimation, we also use it in rank-minimization interpolation as a cost-free initial solution to the optimization problem. As statics estimation and interpolation operate in the midpoint-offset domain, we avoid the cost of transformations back and forth from the source-receiver to the midpoint-offset transform domain. Consequently, our reconstruction, which indicates its potential on synthetic and field data, also is computationally efficient.

Seismic angle gathers and spectral seismic attributes offer complementary insights to improve understanding of complex subsurface characteristics. However, the labor-intensive process of subsurface characterization, data annotation, limited labeled data, and subsurface complexity make it difficult to leverage these insights via supervised learning approaches.

To overcome such challenges and benefit from the strength of spectral seismic attributes, this study introduces a novel hierarchical Self-Organizing Map (SOM) framework to integrate spectral seismic attributes like scalograms and spectrograms (joint time-frequency analyses) extracted from angle gathers.

In our current research, firstly, we trained individual SOMs, as an unsupervised pattern recognition algorithm on reflectivity images, angle-gathers, and the spectral seismic attributes extracted from angle-dependent data. Secondly, we deploy a hierarchical SOM network to combine and analyze all these datasets. Thirdly, we evaluate the hierarchical approach and standalone analyses of clustering quality and information content using the binary boundary maps and the performance metrics. Our findings indicated that, the scalogram-based hierarchical SOM, containing information of different angles, achieves the lowest Quantization Error and Davis-Bouldin Index, indicating optimal feature representation and well-separated clusters. The findings stress the potential of hierarchical networks and joint time-frequency analyses from angle gathers for robust seismic interpretation workflows.

With a focus on geo-imaging applications for the energy transition, we are looking for affordable, but still accurate seismic imaging methodologies. One of those recently developed methods is Joint Migration Inversion, which involves the joint estimation of the seismic reflectivity image and the background propagation velocity model. This method operates in the frequency domain and is based on recursive wavefield propagation, while including all scattering and transmission effects. The involved full wavefield modeling engine is the most time-consuming part of the JMI process, so accelerating this has direct impact on the overall costs. One option is making use of the fact that the modeling can be done independently per frequency component, such that we can model the data for a subset of these frequencies and use interpolation to obtain the data at missing frequencies. This papers studies the use of a neural network (NN) approach for this interpolation process. We investigate the accuracy of the interpolation process under different sub-sampling ratios and using regular or irregular subsampling. The counter-intuitive result is that regular subsampling gives slightly better results. Moreover, we demonstrate that we can go down to 66% missing frequencies with the currently used NN based on the cGAN approach.

The development of many offshore wind-farms around the world is one way to reduce CO2 emissions in the atmosphere. Each new offshore wind-energy site needs to be characterized to find the optimal locations for the monopiles or other wind-turbine infrastructure. For a global investigation this characterization is usually done with seismic surveys, in order to get a good view on the first 100m of the subsurface. The main seismic parameter that can be well defined from acoustic reflection measurements is the P-wave velocity of the subsurface layers, while the subsurface strength is much better coupled to the shear wave velocity or the related shear-wave modulus. One way to obtain such information from marine reflection seismic data is via AVO information embedded in PP reflections. However, the typical low values for the S-wave velocities in the first 10s of meters below the sea-bottom make it not a favorable parameter to estimate. Therefore, based on synthetic data modeling, the feasibility of estimated S-wave velocity information from marine reflection data is investigated. Different scenarios with varying maximum offset are investigated, where the conclusion is that the ability to retrieve accurate S-wave velocity information is limited, but can be improved by acquiring more offsets

Seismic data acquired on land face multiple challenges due to the near-surface complexity. One of the challenges is the weathering layers’ influence due to the low velocity and rapidly varying nature of these layers. To overcome that, dense source-receiver sampling can be used to characterize the near-surface. However, that increases the acquisition costs, which makes it less attractive for large-scale seismic acquisition. An alternative approach is to use compressive sensing to acquire the data. Although it enables economical data acquisition, compressive sensing requires data reconstruction, which is another challenge in the presence of complex weathering layers. To overcome that, we propose a joint reconstruction and near-surface correction algorithm using a model-independent low-rank-based approach. We apply the method to synthetic and real data, which shows superior results compared with the conventional approach of near-surface correction followed by data reconstruction.

Fluvial reservoirs are a major target for geothermal energy production. Interpreting the 3D reservoir architectures from 2D seismic datasets, which usually are acquired for geothermal systems, is difficult. In particular, small-scale geological factors like sandbody connectivity are challenging to resolve. This study addresses these issues through a novel workflow that incorporates 3D geological and 2D seismic modelling methods to assess the seismic responses of stratigraphic attributes in fluvial geothermal reservoirs where data availability is low.

Two synthetic fluvial reservoir scenarios were built, ranging from a single channelised deposit to a geologically more plausible model ensemble of fluvial deposits, which represents the reservoir heterogeneities that could be present at the geothermal doublet at Delft University of Technology. Acoustic finite-difference modelling was combined with seismic imaging to create 2D depth images. Our results reveal how seismic resolution determines our ability to correctly identify sandbody connectivity and capture inner channel details. Whereas channel bodies can be detected, the best frequency spectra for observing certain geological features remain unclear. These findings emphasise that quantitative multi-scale analysis, advanced imaging techniques, and survey design optimisation are central to improving seismic characterisation of fluvial geothermal systems in future research.

This study examines the applicability of seismic methods for monitoring hydrogen storage and detecting potential leakage in sandstone reservoirs, with a particular focus on amplitude variations in angle-dependent image gathers. Using the FluidFlower benchmark model as a controlled geological framework, two types of sandstone—mildly consolidated and unconsolidated—are considered. Gassmann’s fluid substitution is used to model elastic property changes under different hydrogen saturation and leakage scenarios, and seismic responses are generated using Kennett’s reflectivity method.

The analysis shows that seismic amplitudes are sensitive to both fluid saturation and lithology. In mildly consolidated sandstones, hydrogen injection leads to observable increases in amplitude at reservoir interfaces. In unconsolidated sandstones, elastic contrasts are more pronounced, resulting in stronger and more detectable seismic responses. These findings highlight the need to account for lithological characteristics when designing seismic monitoring strategies for underground hydrogen storage.

High-resolution seismic reflections are essential for imaging and monitoring applications. In seismic land surveys using sources and receivers at the surface, surface waves often dominate, masking the reflections. In this study, we demonstrate the efficacy of a two-step procedure to suppress surface waves in an active-source reflection seismic data set. First, we apply seismic interferometry (SI) by cross-correlation, turning receivers into virtual sources to estimate the dominant surface waves. Then, we perform adaptive subtraction to minimize the difference between the surface waves in the original data and the result of SI. We propose a new approach where the initial suppression results are used for further iterations, followed by adaptive subtraction. This technique aims to enhance the efficacy of data-driven surface-wave suppression through an iterative process. We use a 2-D seismic reflection data set from Scheemda, situated in the Groningen province of the Netherlands, to illustrate the technique’s efficiency. A comparison between the data after recursive interferometric surface-wave suppression and the original data across time and frequency–wavenumber domains shows significant suppression of the surface waves, enhancing visualization of the reflections for subsequent subsurface imaging and monitoring studies.

Reflection waveform inversion (RWI) is a technique that uses pure reflection data to estimate subsurface background velocity, relying on evolving seismic images. Conventional RWI operates in a cyclic workflow, with two key components in each cycle—migration and reflection tomography. Conventional RWI may result in suboptimal background velocity estimation, partly due to limited or unresolved resolution within each component in each cycle. While gradient pre-conditioning with the reciprocal of Hessian information helps resolve this issue in both components of RWI, it becomes impractical for a large number of model parameters. One-way reflection waveform inversion (ORWI) is a reflection waveform inversion technique in which the forward modelling scheme operates in one direction (downward and then upward) via virtual parallel depth levels within the medium. Leveraging the ORWI framework, we decompose and reduce the linear Hessian operator (also known as the approximate Hessian or Gauss–Newton Hessian) into multiple smaller suboperators. In particular, the diagonal blocks of the monofrequency approximate Hessian operators, each corresponding to a single depth level within the medium, are extracted and inverted to pre-condition the corresponding monofrequency gradients in both the migration and reflection tomography components of ORWI. This depth-dependent gradient pre-conditioning transforms standard ORWI into a high-resolution, yet computationally feasible version aimed at addressing suboptimal velocity estimation, referred to as high-resolution ORWI. The effectiveness of the proposed approach is demonstrated through successful applications to synthetic data examples.

Reducing the uncertainty of reservoir characterization requires to better identify the small-scale structures of the subsurface from the available data. Studying the seismic response of meter-scale, stratigraphic heterogeneities typically relies on the generation of reservoir models based on outcrop examples and their forward seismic modelling. To bridge geological information and seismic modelling, these methods allocate values of acoustic properties, such as mass-density and P-wave velocity, according to discretized properties like layer-type lithology or facies units. This strategy matches the current workflow in seismic data inversion in industry, where modelling workflows are based on lithofacies distributions. However, from stratigraphic modelling, we know that meter-scale heterogeneities occur within certain facies and lithologies. Here, we evaluate the difference on the seismic response between allocating acoustic properties in a grain size–based, semi-continuous manner versus discretized manners based on lithology and facies classifications. To do so, we generate a reference geological simulation that we populate with acoustic properties, mass-density and P-wave velocity, using three different strategies: (1) based on grain size distribution; (2) based on facies distribution; and (3) based on lithology. The method we propose includes the generation of realistic geological simulations based on stratigraphic modelling and the transformation of its output into acoustic properties, honouring the intra-lithology and intra-facies, small-scale structures. We, then, generate seismic data by applying a forward seismic modelling workflow. The synthetic data show that the grain size–based simulation allows the identification of small-scale, stratigraphic heterogeneities, such as beds with strong density and velocity contrasts. These stratigraphic structures are smoothened or may completely disappear in the facies and lithology discretized simulations and, therefore, are not (well) represented in the synthetic seismic data. Recognizing meter-scale, stratigraphic heterogeneities is relevant for the characterization of the fluid flow in the reservoir. However, current discrete and lithology-based strategies in seismic inversion are not able to resolve such heterogeneities because real subsurface properties are not discrete properties but continuous, unless there are stratigraphic discontinuities such as erosional surfaces or faults. This research works towards a better understanding of the relationship between changes in these continuous properties and the observed seismic data by introducing greater complexity into the discretized geological simulations. Here, we use synthetic seismic images with the goal of eventually aiding in fine-tuning seismic inversion methodologies applied to real seismic data. One pathway is to foster the development of inversion approaches that can leverage stratigraphic modelling to get stronger geological priors and replace the standard but inadequate multi-Gaussian prior.

Seismic wave attenuation, quantified by the quality factor (Q), leads to energy loss and waveform distortion, significantly degrading seismic data quality and resolution. Accurate Q estimation is essential for understanding subsurface properties, particularly in applications such as carbon capture and storage (CCS) and near-surface studies, where attenuation effects are pronounced due to the presence of fluids, gases, or loose soil. Traditional Q tomography methods predominantly rely on spectral-ratio or centroid-frequency shift approaches to account for attenuation effects. However, these methods often face significant limitations, including oversimplified wave propagation assumptions, poor localization in heterogeneous media, and a tendency to produce smeared results, ultimately reducing resolution and accuracy.

To address these challenges, we introduce a novel Q-estimation approach that integrates full-waveform matching for accurate attenuation-effect estimation and compensation during the migration process. The Full Wavefield Migration method is enhanced by incorporating Q into a one-way modeling operator, utilizing full-waveform matching for precise Q estimation, and applying a Random Forest regression constraint to mitigate cross-talk between Q and reflectivity. This approach enables robust and localized Q estimations. Numerical examples demonstrate its effectiveness in accurately retrieving both reflectivity and attenuation models, thereby improving imaging resolution in complex subsurface environments.

The overburden structures often can distort the responses of the target region in seismic data, especially in land datasets. Ideally, all effects of the overburden and underburden structures should be removed, leaving only the responses of the target region. This can be achieved using the Marchenko method. The Marchenko method is capable of estimating Green's functions between the surface of the Earth and arbitrary locations in the subsurface. These Green's functions can then be used to redatum wavefields to a level in the subsurface. As a result, the Marchenko method enables the isolation of the response of a specific layer or package of layers, free from the influence of the overburden and underburden. In this study, we apply the Marchenko-based isolation technique to land S-wave seismic data acquired in the Groningen province, the Netherlands. We apply the technique for combined removal of the overburden and underburden, which leaves the isolated response of the target region, which is selected between 30 and 270 m depth. Our results indicate that this approach enhances the resolution of reflection data. These enhanced reflections can be utilised for imaging and monitoring applications.

As seismic imaging moves towards the imaging of more complex media, properly modelling elastic effects in the subsurface is becoming of increasing interest. In this context, elastic wave conversion, where acoustic, pressure (P-) waves are converted into elastic, shear (S-) waves, is of great importance. Accounting for these wave conversions, in the framework of forward and inverse modelling of elastic waves, is crucial to creating accurate images of the subsurface in complex media. The underlying mechanism of wave conversion is well understood and described by the Zoeppritz equations. However, as these equations are highly nonlinear, approximations are commonly used. The most well-known of these approximations is Shuey’s approximation. However, this approximation only holds for small angles and small contrasts, making it insufficient for realistic forward and inverse modelling scenarios, where angles and contrasts may be large. In this paper we present a novel set of approximations, based on Taylor expansions of the Zoeppritz equations, which we name the extended Shuey approximations. We examine the quality of these approximations to the Zoeppritz equations and compare them to existing approximations described in literature. We then apply these extended Shuey approximations to the elastic full-wavefield modelling algorithm for a simple, synthetic, 1.5-D example, where we show that we can accurately model the P- and S-wavefields in a forward modelling case. Finally, we apply our approximations to the elastic full-wavefield migration algorithm for a simple, synthetic, 1.5-D example, where we show that we can recover an accurate image in an inverse modelling case.

In recent years, deep learning (DL) has emerged as a promising alternative approach for various seismic processing tasks, including primary estimation (or multiple elimination), a crucial step for accurate subsurface imaging. In geophysics, DL methods are commonly based on supervised learning from large amounts of high-quality labelled data. Instead of relying on traditional supervised learning, in the context of free-surface multiple elimination, we propose a method in which the DL model learns to effectively parameterize the free-surface multiple-free wavefield from the full wavefield by incorporating the underlying physics into the loss computation. This, in turn, yields high-quality estimates without ever being shown any ‘ground truth’ data. Currently, the network reparameterization is performed independently for each dataset. We demonstrate its effectiveness through tests on both synthetic and field data. We employ industry-standard Surface-Related Multiple Elimination (SRME) using, respectively, global least-squares adaptive subtraction and local least-squares adaptive subtraction as benchmarks. The comparison shows that the proposed method outperforms the benchmarks in estimation accuracy, achieving the most complete primary estimation and the least multiple energy leakage, but at the cost of a higher computational burden.

Seismic Imaging (SI) survey design for onshore applications faces challenges such as accessibility and poor data quality due to unexpected (near-)surface conditions. In this paper, we explore the correlation between the surface conditions provided by Land-Cover (LC) maps generated using Remote Sensing (RS) data and different settings of Seismic Processing (SP) parameters. The study involves a 2D seismic line related to geothermal exploration from the Netherlands.

The accuracy of a model obtained by full-waveform inversion can be estimated by analysing the sensitivity of the data to perturbations of the model parameters in selected subsurface points. Each perturbation requires the computation of the seismic response in the form of Born scattering data for a typically very large number of shots, making the method time consuming. The computational cost can be significantly reduced by considering the point where the subsurface parameters are perturbed as a Born scatterer. Instead of modelling each shot separately, reciprocity relations provide the Green functions from the sources to the scatterer in terms of Green’s functions from the scatterer to the sources. In this way, the Born scattering data from a single point in the isotropic elastic case for a marine acquisition with pressure sources and receivers can be expressed in terms of the Green functions for force and moment tensor sources located at the scatterer and only a small number of forward runs are required. A 2-D example illustrates how the result can be used to determine the hessian and local covariance matrix for the model parameters at the scatterer at the cost of 5 forward simulation.

The data-driven surface-related multiple elimination (SRME)-type approach requires fully sampled sources and receivers during the multidimensional convolution process. Otherwise, the estimated multiples will be aliased. Compared to expensive reconstruction processes before prediction, dealiasing on the estimated multiples from limited sources might provide a potential easier solution in a 2D scenario, where deep learning (DL)-based methods suit well for this highly non-linear problem. Unfortunately, DL-based multiple dealising will not function well for 3D data due to extremely coarse sampling in either source or receiver side. Thus, data interpolation/reconstruction is the only option, though the performance might not be desired. Generalized surface multiple prediction (GSMP) is the most used on-the-fly interpolation approach in 3D. Still, GSMP accuracy heavily relies on the existing traces. When fed with coarsely sampled recorded data only, GSMP tends to generate multiples with low amplitude and distorted phase, especially for small offsets. We propose a U-Net framework to repair GSMP estimated multiples such that the amplitude loss and distorted phase can be restored. In this way, the strong non-linear mapping power from DL can help repair the GSMP estimated multiples.