Geophysical Prospecting - Volume 73, Issue 9, 2025

Accurately and effectively handling undulating interfaces, including both free surfaces and internal interfaces, remains a key challenge in seismic exploration. Traditional scalar wave equations typically neglect the influence of internal interfaces on wave propagation. To address this, the wave equation in layered media (WEILM) is established by introducing the Dirac delta function. For undulating interfaces, existing methods are mostly grid‐based, and often require additional grid processing to achieve accurate description, which increases computational cost. Therefore, by introducing the meshless method combined with free surfaces boundary conditions, this article proposes the method for dealing with undulating surfaces on the basis of the radial‐basis‐function‐generated finite difference (RBF‐FD) method. Theoretical analysis indicates that the stability of the proposed method is affected by the sum of the stencil node weights. Numerical experiments show that, compared with the scalar wave equation, the interfaces term introduced in WEILM effectively adjusts waveform amplitudes. Moreover, relative to the classical Lax–Wendroff correction (LWC) method, our approach can avoid spurious diffraction waves caused by grid discretization when handling undulating surfaces. By applying RBF‐FD to solve WEILM in conjunction with the method for dealing with undulating surfaces, complex seismic wavefield, including undulating surfaces, can be simulated with higher accuracy.

The limited understanding of rock physical properties in marine shale from the Ordovician Wulalike Formation along the western margin of the Ordos Basin has hindered comprehensive evaluations of shale gas reservoirs. This study systematically investigates the variation patterns and controlling factors of seismic elastic properties in Wulalike Formation marine shale samples through petrological and petrophysical tests, while discussing the distribution characteristics of reservoir ‘sweet spots’. Results indicate that the petrological characteristics of Wulalike Formation marine shale are influenced by tectono‐sedimentary differentiation. The lithology transitions from calcareous shale in upper slope environment to mixed shale in slope depression, and finally to siliceous shale in open marine shelf environment. Both organic matter abundance and porosity of the shale samples progressively increase with depositional environments and lithological transitions. Simultaneously, the rock stress skeleton evolves from carbonate particle dominance to clay and quartz particle dominance. Variations in rock microstructural characteristics among different lithological types of samples are the primary factor influencing seismic elastic properties. In petrophysical crossplots (impedance vs. porosity, Poisson's ratio vs. P‐wave impedance and λρ vs. μρ), the shale samples exhibit partitioned distributions on the basis of their composition and lithology. The ‘sweet spots’ reservoirs are predominantly composed of siliceous shale, characterized by high total organic carbon (TOC), porosity and low Poisson's ratio and λρ characteristics. On the basis of the petrophysical analysis, reservoirs are categorized into three grades. Laterally, reservoir classification transitions from Grade ‘III’ to ‘I’ with changing depositional environments. Shale gas reservoirs in open marine shelf and slope depression environments (e.g., the lower part of Well ZP1) meet or exceed Grade ‘II’ standards, indicating high‐quality reservoir potential.

Seismic fluid discrimination plays a critical role in sweet spot detection, reservoir characterization, reserve evaluation and well placement. Tight sandstone reservoirs are typically characterized by low porosity, poor pore connectivity, complex pore types, non‐uniform gas–water distribution and strong heterogeneity, which often lead to inaccurate fluid discrimination. In this study, we develop a double‐porosity equivalent medium model for tight sandstone reservoirs using the Keys–Xu model combined with Gassmann's equation. We systematically investigate the effects of pore structure, porosity and water saturation on elastic responses. On the basis of this model, a rock physics template (RPT) is constructed using the P‐wave modulus and the P‐ to S‐wave modulus ratio. Polynomial fitting is then applied to derive mathematical expressions for both water‐ and gas‐saturated trendlines. On the basis of these trendlines, an RPT‐based fluid indicator is defined to quantify deviations from the gas‐saturated sandstone trendline. We further apply the proposed fluid indicator to a tight gas sandstone reservoir in the central Sichuan Basin, Southwest China. The strong agreement between the extracted fluid indicator and well log‐based water saturation interpretation demonstrates that this method significantly improves the accuracy of fluid content quantification compared with traditional semi‐quantitative RPT‐based approaches. Application to seismic data further shows that our method yields a reasonable estimation of gas distribution in tight sandstone reservoirs, confirming its reliability and practical applicability for fluid characterization. This approach offers promising potential for quantifying fluid content in deep‐buried tight reservoirs.

Pre‐ and post‐stack seismic inversion is the primary approach for converting collected seismic data into geophysical property models particularly velocity for subsurface interpretation and reservoir characterization. The traditional workflows often start from an initial property model and iteratively revise it by minimizing the misfit between the acquired real seismic dataset and the synthetic one derived from the updated property models, here denoted as the soft constraint in seismic space. However, it heavily relies on human supervision in building a good initial model and monitoring the misfit optimization process. On the contrary, most of the recent deep learning‐based workflows target at non‐linearly mapping seismic patterns to properties measured at available well only, here denoted as the hard constraint in property space. Correspondingly, its accuracy greatly depends on the availability of sufficient training wells; otherwise, overfitting occurs causing the machine prediction not meeting the soft constraint throughout the target seismic survey. To resolve these limitations, this study presents a practical workflow that enables rock property inversion from partial‐stack seismic via two physics‐guided convolutional neural networks (CNNs), with the first one embedding the approximated AVO gradient to build initial property models that satisfy the hard constraint and the second one embedding the reflection coefficients to refine the models by enforcing the soft constraint. Between both CNNs is the use of well‐established relevant physics to generate pseudo property–reflectivity–seismic pairs for training the second CNN. Its added values are validated through applications to the Volve survey in North Sea and the Exmouth survey in Western Australia. The produced property volumes not only are observed of high lateral consistency and vertical resolution but also derive synthetic seismic data that are closely correlated with actual seismic data.

Hydrothermal alteration plays a crucial role in the precipitation of gold and other metals, particularly within orogenic gold deposits hosted in mafic and ultramafic rocks. This alteration significantly modifies the rock matrix, leading to changes in its petrophysical properties. In this study, we focus on two key processes: quartz‐carbonate vein formation and sulfidation, both of which have distinct effects on geophysical measurements. The investigation centres on the Piché Group within the Augmitto–Bouzan sector of Rouyn Property, a primary target for gold exploration. Quartz‐carbonate vein formation, characterized by the direct precipitation of resistive minerals such as quartz and carbonates, has a pronounced effect on resistivity logs, leading to increased resistivity values. These changes are also evident in sonic logs, where P‐wave and S‐wave velocities increase due to the presence of these minerals. Sulfidation, in contrast, reflects metasomatic alteration of the host rock and is primarily captured through induced polarization (IP) and spontaneous potential (SP) logs. The crystallization of sulfide minerals, such as pyrite and arsenopyrite, not only leads to increased IP and decreased SP values but also results in a high degree of variability in IP values, reflecting the heterogeneous distribution of sulfides in the host rocks. These findings are further supported by micro‐XRF data, confirming the presence and distribution of sulfide minerals in key alteration zones. Our results suggest a relationship between these petrophysical signatures and gold endowment, with increased resistivity (presence of quartz‐carbonate veins) and elevated IP values (sulfidation) correlating with a higher probability for gold concentrations. By focusing on the distinct effects of quartz‐carbonate veins and the sulfidation process, this study provides valuable insights into the identification of alteration zones and their potential for mineralization. These results contribute to a deeper understanding of alteration within the Piché Group and offer a framework for more targeted exploration efforts in similar geological environments.

The transformative impact of deep‐learning architectures on machine learning has been substantial. Recently, a wide range of studies have successfully applied these methods to seismic facies segmentation using well‐established public datasets, such as F3 and SEAM AI. However, many of these works lack detailed descriptions of their methodologies and implementation details, including dataset partitioning, hyperparameter settings and other critical aspects. The lack of reproducibility information makes fair comparison between studies quite difficult, as methodological details can heavily affect the results obtained. In this work, we discuss this problem and present a fair comparison between five state‐of‐the‐art models commonly used in the literature: DeepLab V3, DeepLab V3+, Segmenter, SegFormer and SETR. We found that the SETR model has promising performance on both the F3 and SEAM AI datasets and convolutional neural network models offer a higher performance to parameter count ratio compared to the transformer models.

Interpretation of magnetotelluric data acquired near coastlines is challenging due to the distortion of electromagnetic fields caused by the sea effect. Specifically, it is essential to accurately simulate the propagation of electromagnetic fields around the land–sea and air–sea boundary. To properly address the sea effect from the shallow sea, we present a hybrid edge‐based finite element method that combines prismatic and tetrahedral elements with vector shape functions. The 3D mesh incorporates prismatic elements to stably obtain enough vertical resolution for the land–sea and air–sea boundary, whereas tetrahedral elements are used for thicker subsurface volumes. The proposed method is demonstrated using synthetic data from a simple 1D/2D model for which analytic/pseudo‐analytic solutions can be obtained and used as reference data. These examples show that the proposed method achieves computational efficiency while maintaining accuracy across domains with thin geometries.

Methane is a primary component of natural gas that is the cleanest fossil fuel to support the sustainable development of our society. Electrical survey methods are frequently employed for the exploration of methane that is majorly existing in hydrocarbon systems in the subsurface earth. However, although the quantitative interpretation of electrical survey data relies on the knowledge about the electrical properties of methane bearing rocks, it remains poorly understood about how dissolved methane affects the electrical resistivity of brine and brine‐saturated sandstones. We bridge this knowledge gap by measuring the electrical resistivity of brine and brine‐saturated artificial sandstone samples with varying brine salinity, pore pressure and temperature, as a function of dissolved methane. We find that the dissolution of methane improves the electrical resistivity of brine, and the improvement can be best fitted by a regression equation comprising both a constant and an exponential part. We also find that the improvement in the brine resistivity increases with reducing brine salinity, pore pressure and temperature, where reducing brine salinity, pore pressure and temperature also improve the brine resistivity with no dissolve methane. Experiment on the rock samples shows that the resistivity of the artificial sandstones with dissolved methane behaves in a similar way to the brine resistivity, in terms of its dependence on the brine salinity, pore pressure and temperature. Further analyses demonstrate that the dissolution of methane in brine does not affect the cementation exponent of the rocks, and therefore the rock resistivity with dissolved methane can be predicted on basis of its constant cementation exponent and the brine resistivity with dissolved methane. The results not only reveal the effects of dissolved methane on the electrical resistivity of brine and sandstones at varying conditions but also pave the way to the interpretation of electrical survey data for the better quantification of dissolved methane in the hydrocarbon systems.

Stochastic reservoir characterisation relies on the careful integration of geological modelling and geophysical data, enabling prediction and uncertainty quantification for reservoir decision‐making. In this paper, we address some of the computational challenges associated with Bayesian reservoir characterisation, focusing on key obstacles: demanding geophysical forward modelling, high dimensionality of spatial variables and effective posterior sampling of reservoir variables given geophysical data and well information. Leveraging a pseudo‐Bayesian approach, we replace the intricate forward model for seismic amplitude‐versus‐offset data with a computationally efficient multivariate adaptive regression splines method, resulting in a 34‐times acceleration in computations. For handling high‐dimensional variables modelled by Gaussian random fields, we employ a fast Fourier transform technique. We use a preconditioned Crank–Nicolson method for efficient Markov chain Monte Carlo sampling from the posterior of the reservoir variables. The approach is motivated by challenging reservoir conditions at the Alvheim field in the North Sea, where we demonstrate our approach for Bayesian posterior sampling of oil and gas saturation and clay content conditional on seismic amplitude data and well information. We compare our results with those obtained from an approximate ensemble‐based Kalman method for posterior sampling.

Seismic data interpolation using convolutional neural networks (CNNs) suffers from accuracy limitations due to the inter‐band interference across different frequency bands, which negatively affects subsequent inversion and interpretation. To address this limitation, we propose a multi‐band strategy that first decomposes the seismic data into multiple sub‐bands through frequency filtering. Independent CNN models are then used to process each specific frequency band to isolate spectral interference. We focus on regularly missing shots interpolation, assuming that dense receiver arrays are available during seismic acquisition with sparse shots. As for the training data preparation, the spatial reciprocity of Green's function is considered, which guarantees the similarity between common shot gathers (CSGs) and common receiver gathers (CRGs). The available dense CSGs are used to train networks using the multi‐band‐assisted training strategy. The resulting optimized independent models are then employed to reconstruct missing shots in sparse CRGs for each frequency band separately. Interpolated multi‐band data are finally fused by summation to obtain the full‐band result. Numerical experiments on synthetic and field data demonstrate that the proposed multi‐band‐assisted training strategy provides superior interpolation accuracy compared to traditional full‐band training, particularly in mitigating cross‐band interference.

Missing seismic traces from data acquisition limits often significantly degrade data quality. This study presents an unsupervised method using implicit neural representation (INR), specifically sinusoidal representation network (SIREN), to enhance seismic data quality from a single shot gather. Notably, the unsupervised framework trains the SIREN by optimizing it on the observed traces in the single‐gather data. The network learns a continuous function, enabling the reconstruction of missing data at any spatio‐temporal coordinate. This algorithm directly addresses both missing trace interpolation and the enhancement of sparsely sampled data resolution. Key network design choices, such as exponential frequency scaling and dense skip connections, are shown to enhance reconstruction accuracy by mitigating spectral bias and incorporating multi‐scale features. Furthermore, our analysis of different coordinate handling strategies identifies a key trade‐off on the geometry setting. Reframing interpolation as a super‐resolution task enables the successful reconstruction of up to 75% regularly missing traces and can maintain continuity across large gaps of up to 10 traces. However, this method proves geometrically inaccurate for irregular missing data, as it discards true physical coordinates, leading to incorrect solutions. In contrast, strategies that maintain physical coordinates show significantly degraded performance when faced with such large‐scale data gap. The proposed framework successfully interpolated multichannel seismic data and enhanced sparse ocean bottom cable (OBC) data resolution. Although challenges remain for large irregular gaps and computational efficiency, this work establishes SIREN as a promising unsupervised tool for single‐gather seismic interpolation and sparse data resolution enhancement without requiring external training data.

Staggered‐grid finite‐difference (FD) schemes are widely used in numerical simulation of seismic wave propagation. The traditional explicit staggered‐grid scheme adopts the second‐order temporal and explicit high‐order spatial FD, so it easily suffers from significant temporal dispersion and the spatial operator–length saturation effect. The recently developed implicit time–space‐domain high‐order staggered‐grid scheme for 2D acoustic wave simulation overcomes those two weaknesses effectively, yielding high‐order accuracies in both time and space and thus better suppressing the numerical dispersion. However, the involved FD coefficients are generally determined by Taylor‐series expansion (TE) or the least‐squares (LS) method and still cannot effectively control spatial dispersion at the large wavenumber range. Adopting the same FD stencil, we alternatively determine the FD coefficients using a combination of TE and the Remez optimization algorithm. The temporal accuracy–related coefficients are determined by the TE of the time–space‐domain dispersion relation, whereas the implicit spatial FD coefficients are calculated by using the Remez exchange optimization algorithm to effectively extend the effective wavenumber range while achieving a high‐order temporal accuracy and consequently enhance the overall modelling accuracy. The newly optimized scheme is then extended into a 3D case. Dispersion and numerical analyses validate that the proposed new schemes better suppress the spatial dispersion while maintaining the high‐order temporal accuracy and outperform the existing TE‐ and LS‐based FD schemes. To further improve the modelling efficiency, the variable–operator–length strategy is combined. Numerical examples of the 2D and 3D complicated models validate the effectiveness of the combined scheme in reducing the operator length and consequently improving the modelling efficiency.

This study describes an integrated landslide monitoring program in the landslide‐prone Tizzano Val Parma region (Italy) using traditional and innovative geoelectrical techniques, namely, electrical resistivity tomography (ERT), induced polarization (IP) and self‐potential (SP) methods. Both the conventional fixed‐base and the unconventional sparse gradient array configurations were adopted. The use of analytic signal amplitude (ASA) technique enabled for a better recognition of primary SP anomaly sources, for the sparse gradient arrays, providing useful insights in delineating areas of interest. The region faces recurrent landslides due to geological and geomorphological factors, leading to high hydrogeological instabilities and environmental risks. Borehole stratigraphy reveals a complex lithology of sandstones, clayey marl, and coarse materials. The ERT–IP survey provides insights into various landslide types, identifying distinct domains including complex quiescent, active and undetermined landslides. Active fault evolution is observed, indicating potential risk zones. Sparse gradient SP monitoring captures short‐term electrokinetic anomalies and stable long‐term variations between 2022 and 2023. Both fixed‐base and sparse gradient SP monitoring highlight displacement anomalies towards the valley, suggesting potential landslide movements. Interpretation of SP maps allowed to identify preferential water flow directions which denote probable accentuating risks during intense rainfall events. This study emphasizes the significance of integrated geoelectrical monitoring for early landslide detection. Non‐conventional and conventional SP arrays provide insights into anomaly repeatability and stability. Comparison of ERT–IP results with borehole information enables the extrapolation of geological characteristics, yielding a holistic understanding of subsurface structures and potential risk zones. Integration of these methodologies contributes to effective landslide management, underscoring the dynamic nature of landslide‐prone regions and the necessity for ongoing risk assessment and monitoring.

Deep learning has been widely used to invert geophysical properties due to the availability of training data and an increased computing power. In particular, Bayesian deep learning is commonly applied to estimate the uncertainty of rock properties, which is essential for risk management and decision‐making. However, the selection of appropriate prior parameters, such as the standard deviation in the Gaussian prior distribution placed on neural parameters including neural weights and biases, is crucial for training Bayesian neural networks (BNN), as it significantly impacts the prediction performance of the trained models. In this research, we introduce a hierarchical structure to the BNN, and the mean and standard deviation in the Gaussian prior placed on neural parameters are randomly drawn from hyper‐priors, thus excluding the preliminary tuning runs with trial values. Compared to traditional BNN, a consistent prediction accuracy is achieved with an estimate of aleatoric and epistemic uncertainties when using the hierarchical Bayesian networks with different hyper‐priors, thereby making them more robust against the choice of prior parameters. In addition, we apply the statistical sampling technique to reduce the overall size of the training data, which can proportionally decrease the training time of the deep learning models when large amounts of training data are available.

Three‐dimensional (3D) diffraction processing aims at superresolution by imaging small‐scale geological features of the subsurface localized as points and space curves. In analogy to the (anti‐) stationary phase filtering, we separate images of points from images of lines by weighting the Kirchhoff migration. In addition to the deviation from the specularity and Snell's law, the new summation weights verify the conformity of seismic traces to Keller's law of edge diffraction. In addition to that, the configuration of the reflectors determines the diffraction phase reversal pattern specific to isolated lines, edges and wedges. To counteract the summation of the opposite phases in 3D, we provide extra alternating factors for edge and wedge diffraction. All these weights require local orientation of diffractors and reflectors, which we simultaneously retrieve from the full‐wave image by a modification of the slant‐stack search. Synthetic examples show the benefits of the proposed techniques.

Vertically transverse isotropic (VTI) acoustic wave equations are widely used to simulate wave propagation in VTI media. A commonly used acoustic VTI wave equation can be derived by setting the vertical shear‐wave velocity to zero in the elastic VTI wave equation. However, the resulting acoustic VTI wave equation has a non‐symmetric propagation operator, which leads to the operator of the adjoint equation in full‐waveform inversion (FWI) being different from that of the forward equation. Consequently, two separate sets of code are required for simulating the forward and adjoint wavefields. To simplify code implementation, we propose a symmetric‐form acoustic VTI equation. This new formulation allows both the forward and adjoint equations in FWI to share the same operator, enabling a unified code for both the forward and adjoint wavefield simulation and streamlined implementation. In addition, although both the symmetric and non‐symmetric formulations yield the same gradient, the adjoint wavefield from the non‐symmetric equation shows weaker amplitudes in deeper regions compared to that from the symmetric equation. As a result, FWI using the non‐symmetric formulation may suffer from insufficient compensation when employing a spatial preconditioner based on an approximated diagonal pseudo‐Hessian, leading to slower convergence. Numerical examples using the Marmousi2 and BP anisotropic models, as well as an ocean‐bottom cable (OBC) field data set, demonstrate that the proposed symmetric‐form acoustic VTI FWI achieves better inversion results and faster convergence than its non‐symmetric counterpart.