Understanding fault-zone permeability is crucial in model-based assessment of fluid migration, earthquake nucleation, and hydrothermal or hydrocarbon systems. Vertical seismic profiling (VSP) often captures Stoneley (tube) waves generated by fluid-formation coupling in and around a borehole. Tube waves offer valuable constrains to local hydraulic properties. Full simulation of the generation of tube waves using Biot's poroelastic equations is very important, but computationally demanding due to the multiscale nature of the problem, involving fine-scale borehole geometry and long-wavelength seismic wave propagation in the layered media. We develop a semi-analytical approach that can predict borehole pressure response of a normally incident plane P wave in layered poroelastic media, including irregularities in the borehole radius. The model accounts for three key mechanisms for tube-wave generation: (a) due to elastic impedance contrasts, (b) due to fluid infiltration from poroelastic layers, and (c) due to borehole-radius changes. Using a propagator-matrix formulation under low-frequency assumptions, we derive closed-form expressions for the tube-wave amplitudes and validate them using finite-difference poroelastic simulations. The results show that elastic boundaries produce tube waves with opposite polarities, while a thin porous layer and a thin elastic layer generate asymmetric responses with notably different frequency spectra. Our approach improves upon previous effective-source models by accounting for the tube-wave velocity contrasts and ensuring the consistency with the poroelastic theory. This efficient modeling framework enables clearer interpretation of VSP data in fault zones, providing insights that aid in quantitative estimation of the local hydraulic properties.

Plain Language Summary
Underground rock formations, especially in fault zones, can store and move fluids in ways that affect the occurrence of earthquakes and the extraction of subsurface energy resources. Understanding how easily a fluid flows through these rocks, known as permeability, is vital to earthquake science, geothermal energy, and groundwater management. One established approach to investigate permeability uses seismic waves recorded in water-filled borehole, where characteristic signals called tube waves are sensitive to the surrounding rock and fluid properties. In this research, we develop a new and efficient modeling framework to better understand how tube waves are generated when seismic waves interact with a borehole. Our method uses a simplified mathematical model to calculate how different factors contribute to the tube-wave signals. We study how three main causes—changes in the rock layers, fluid movement in the porous zones, and irregular borehole shapes—affect the tube-wave characteristics. Although the computation is extremely fast, the derived results match closely with those obtained through very complex and computationally expensive simulations. The new approach enables more efficient analysis of the field data and improved capacity to monitor underground fluid pathways, especially in active fault zones.

Although seismic methods using S waves can offer high-resolution images of the shallow soil layers, the use of body S-wave tomography for near-surface water monitoring remains underexplored, and the quantitative interpretation of any observed changes in S-wave velocity (VS) in the field conditions is challenging. We conduct a time-lapse S-wave tomography experiment on a field-scale test dike with controlled water levels, allowing for detailed examination of how VS responds to water infiltration. Our results demonstrate that VS decreases progressively, starting from the high-water-side slope and extending across the dike, as the water level rises, with the most significant changes occurring in the sand body and not in the clay cover. The maximum reduction in VS is approximately 40-60 m/s, corresponding to approximately 25%-30% reduction from the initial condition. We use the squared velocity ratio to evaluate the relative contributions of bulk density and shear modulus to VS changes. In the initially unsaturated zone, both factors contribute significantly to the observed VS changes as the zone becomes fully saturated. In fully saturated zones, we assess the changes in the effective stress using the squared VS ratio. Although the low-water side of the dike shows stress changes that are consistent with numerical modeling, the high-water side shows larger stress changes than expected, possibly due to excess pore pressure during the dynamic flow conditions. These findings highlight the potential of body S-wave tomography for high-resolution, near-surface hydrologic monitoring, and provide insights into the complex interactions between physical properties that influence VS changes under varying water levels in field environments.

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.

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.

Seismic interferometry (SI) retrieves the Green function between two receiver locations using their recordings from a boundary of sources. When using sources and receivers only at the surface, the virtual-source gathers retrieved by SI contain pseudo-physical reflections as well as ghost (non-physical) reflections. These ghost reflections are the results of the cross-correlation or auto-correlation (AC) of primary reflections from two different depth levels, and they contain information about the seismic properties of specific layers in the subsurface. We investigated the application of ghost reflections for layer-specific characterization of the shallow subsurface using SI by AC. First, we showed the technique's potential using synthetic data for a subsurface model with a lateral change in velocity, a gradient in depth for velocity, a thickness change and a velocity change of the target layer. Then, we applied the technique to shallow subsurface field data. We also focused on improving the retrieval of ghost reflections by removing the free-surface multiples and muting undesired events in active-source gathers before applying SI. Our results demonstrate that the ghost reflections can be used advantageously to characterize the layer that causes them to appear in the results of SI. Consequently, they can also provide valuable information for imaging and monitoring shallow subsurface structures.

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 are used to redatum wavefields to a deeper level in the subsurface. 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 elimination of the overburden and underburden. Our results indicate that this approach enhances the resolution of reflection data. These enhanced reflections can be utilised for imaging and monitoring applications.

The absence of information on lateral variability in the soil is detrimental to estimating accurately the local site response in the event of an earthquake. To address this problem, the use of densely sampled seismic data together with sparsely distributed but detailed vertical soil profiles obtained from cone penetration tests (CPTs) is advantageous. This study explores the adaptation of suitable machine learning (ML) approaches to derive reliable, site- and depth-specific correlations between seismic shear-wave velocity (Vs) and cone-tip resistance (qc). Such correlation could be successfully established by combining information from seismic CPT surveys with available borehole information for the Groningen region in the Netherlands. It is found that, even over substantial distances, ML-based techniques offer site- and depth-specific correlations between Vs and qc.

Capturing the spatial variability in soil is crucial for ground response analyses in the context of seismic hazard mitigation. The lateral variability in thickness and properties of the different soil layers is one of the main factors that determines the variability of the ground motion spectrum from one location to another. The absence of such lateral variability information in the subsoil in between the locations of Cone Penetration Tests (CPTs) may be compensated by the use of more densely sampled seismic data. In this research we aim to derive a shear-wave velocity field through seismic full-waveform inversion that yields a model resolution approaching that of high-resolution seismic CPT surveys. Following this, a datadriven correlation between geophysical and geotechnical information is attempted through the application of new machine-learning-based approaches.

The Azambuja fault is a NNE trending structure located 50 km north of Lisbon, the capital and most populous city of Portugal. The fault has been considered as a possible source for the historical, large earthquakes. Understanding this fault is a priority in seismic hazard evaluation of this region. The fault has a clear morphological signature. Miocene and Pliocene sediments are tilted eastward and cut by steeply dipping mesoscale fault segments, presenting reverse and normal offsets with a net downthrow to the east. Neotectonic studies indicate a Quaternary slip on the fault of 0.05–0.06 mm/year. However, no direct evidence of the Azambuja fault affecting the Pleistocene or Holocene sediments was found so far. Here, we present the findings from high-resolution seismic reflection studies using both P- and S-waves over the Holocene deposits. The detection of small-throw faulting in ductile sediments is a challenging task. We show that multiple signatures, like perturbations in the reflection hyperbolae visible in shot and CMP gathers, interruptions of reflectors in stacked sections, lateral seismic velocity variations obtained by horizon velocity analysis, all at coincident locations, strongly suggest that the activity of the Azambuja fault has affected the Holocene sediments in the study area. The lateral velocity variations are corroborated by wavepath eikonal traveltime tomography and velocity analysis supported by seismic modeling. By means of 2D viscoelastic modeling, we explain the absence of fault-related diffractions and negligible back-scattered energy from the fault. Using data from nearby boreholes, we find that the 15 ka old alluvium cover has indeed been disturbed by the presence of shallow fault strands. Considering the estimated vertical throws and the empirical relationships between fault length, co-seismic rupture and magnitude, a slip rate of 0.07 mm/y, slightly larger than previously thought, is expected for this fault.

Bulk-density ( ρ) of soil is an important indicator of soil compaction and type. A knowledge of the spatial variability of in situ soil density is important in geotechnical engineering, hydrology and agriculture. Surface geophysical methods have so far shown limited success in providing an accurate and high-resolution image of 3-D soil-density distribution. In this pursuit, 3-D seismic full-waveform inversion (FWI) is promising, provided the robustness and accuracy of density inversion via this approach can be established in the near-surface scale. Ho wever , simultaneous reconstruction of ρand seismic wave velocities through multiparameter FWI remains a challenging task. Near-surface seismic data are commonly dominated by dispersive surface waves whose velocities are controlled by the value and distribution of shear-wave velocity ( V S ). One major difficulty in estimating reliab ly ρfrom near -surface seismic data is due to the relati vel y low sensiti vity of the seismic w av efield to ρcompared to seismic v elocities. Additionally, the accuracy of the estimated ρdecreases due to error in V S -an issue known as parameter coupling. Parameter coupling makes it difficult to estimate accurately ρwithin the framework of conventional gradient-based FWI. More sophisticated optimization approaches (e.g. truncated Newton) can reduce the effect of parameter coupling, but these approaches are commonl y not af fordab le in near -surface applications due to heavy computational burden. In this research, w e ha v e inv estigated how choosing correctly the force direction of the seismic source can contribute to a higher accuracy of ρestimates through 3-D FWI. Using scattered wavefields, the Hessian, and inversion tests, an in-depth and systematic investigation of data sets corresponding to different force directions has been carried out. A comparison of the scattered wavefields due to a point-localized ρperturbation for different force directions shows the robustness of the horizontal-force data set to noise compared to the vertical-force data set. Fur ther more, for a point-scatterer model, an analysis of the gradients of the misfit function using the Hessian shows that utilizing a horizontal-force source enables one to reconstruct the high-resolution gradient with relati vel y small parameter coupling. Finally, inversion tests for two different subsoil models demonstrate that 3-D FWI on a horizontal-force-source seismic data set is capable of providing a more accurate 3-D ρdistribution in soil compared to a vertical-force-source data set. Our results show that the use of a horizontal-force source might allow avoiding computationally demanding, costly optimization approaches in 3-D FWI.

Intersecting faults are often ignored in the geomechanical simulation of induced seismicity. To investigate the effects of fault intersection and the resulting reservoir geometry on induced seismicity, caused, for instance, by gas extraction, we have developed 3D geomechanical models considering two intersecting normal faults and the surrounding horst structure. We simulate the stress field and the dynamic fault reactivation in a uniformly depleted reservoir. We observe that a smaller intersection angle increases the incremental Coulomb stress at the lower reservoir juxtaposition, thus changing the temporal rupture pattern of the seismic event. In our dynamic simulation, the rupture propagates from the main fault to the secondary fault. We conclude that the fault intersection has important effects on the induced seismicity and should be taken into account when evaluating the seismicity risk in a specific region.

Using post-critical reflection data, it is possible to obtain useful information that allows more reliable geological characterization of the subsurface. However, the strong distortion caused by the phase shift in post-critical wavelets makes the use of post-critical reflections rather challenging. For this reason, an approach which is capable of estimating the phase shift of each wavelet of a reflection event in a data-driven manner is desirable. In this vein, in case the frequency spectrum of a wavelet can be correctly estimated, it is possible to estimate the instantaneous phase shift. In this work, we propose an approach which can perform such estimation based on spectral recomposition of seismic data. We design an inversion approach in order to reconstruct the seismic spectrum of the wavelets of a reflection event, which subsequently allows us to estimate the instantaneous phase of each wavelet of the near-surface reflection events without performing prior velocity analysis and/or critical-angle estimation. After finding the instantaneous phase for each wavelet of a reflection event, we show next how one can find the respective phase shifts that can then be corrected.

The structural health of historic quay walls needs to be evaluated well in light of the new conditions they are subjected to. For that, information about their current subsurface structure and condition of their subsurface constructional elements, but also information about the surrounding subsurface structure is crucial. Such information can be supplied by seismic imaging and characterization. We show preliminary results from a high-resolution S-wave survey we performed at a historic quay wall in Overamstel, the Netherlands. We recorded data along four lines – two parallel and two perpendicular to the quay wall. We used a sledge-hammer and a beam as a source and 10-Hz horizontal-component geophones, both oriented in the crossline direction. We show that applying simple processing along the two parallel lines to obtain stacked sections already allowed extracting useful structural information of the subsurface.

Using the critical angle information of a reflection event, it is possible to calculate several essential physical parameters that are key to reliable geological characterization of the subsurface. However, estimation of the critical angle usually requires several steps of seismic processing. For this reason, an approach which is capable of estimating the critical angle directly from the data is of interest. Once the critical angle is estimated, it is possible to estimate further the Poisson's ratio and the seismic velocities. In this work, we propose an approach which can perform this estimation, based on spectral recomposition of seismic data. We design an inversion scheme in order to reconstruct the seismic spectrum of wavelets of a reflection event, which subsequently allows us to estimate the critical angle of near-surface reflection events without performing prior velocity analysis. After finding the critical angle, we show next how to estimate the Poisson's ratio and the compressional- and shear-wave velocities of the medium above the reflector. The approach leads to quite accurate values for Poisson's ratio even for noisy data, in case the number of layers is not large.

Knowing the location and characteristics of shallow subsurface structures like tunnels, cavities, archeological ruins, etc. is of importance for different disciplines and application. To image and/or characterize such objects of interest, different geophysical methods are used. For imaging of a very shallow network of tunnels, the high-resolution seismic method with active sources provide valuable information. We show the results of analysis of an S-wave profile recorded over a network of very shallow tunnels in the Netherlands. The survey used a high-frequency vibratory S-wave source and horizontal particle-velocity geophones, both oriented in the crossline direction, along three lines. We process the reflection data along one of the lines to obtain a stacked section in depth. We also use a method inspired by seismic interferometry to localize a scatterer along the line. We show that both techniques image well the subsurface structures taking into account the 3D ambiguity of processing 2D data.

To investigate the physical processes behind induced seismicities due to, for example, production of hydrocarbons from a reservoir, most of the earlier studies performed geomechanical simulations on a simple reservoir geometry. The effect of fluid depletion is, in general, simulated for such a simple geometry. Neglecting the contribution of realistic 3-D reservoir geometries can lead to a wrong estimation of the incremental stress field. A reliable estimate of the induced stress field is key to producing meaningful simulation results. We perform geomechanical simulations on a simple fault model as well as a more realistic model based on the known geological structures at the earthquake source-region in Zeerijp region, the Netherlands. Our results demonstrate that the angle of the fault intersection affects the incremental stress field, including the effective normal stress, the shear stress, and hence, the Coulomb stress and the SCU value. Our results also show a shift in the rupture pattern and the location of the maximum slip on the fault plane. We conclude that, to properly evaluate the effects of production activities and to simulate precisely the in-situ stress field and the induced seismicity, the incorporation of a realistic reservoir structure in modelling is essential.

Traditional least-squares full-waveform inversion (FWI) suffers from severe local minima problems in case of the presence of strongly dispersive surface waves. Additionally, recorded wavefields are often characterized by amplitude errors due to varying source coupling and incorrect 3D-to-2D geometrical-spreading correction. Thus, least-squares FWI is considered less than suitable for near-surface applications. In this paper, we introduce an amplitude-unbiased coherency measure as a misfit function that can be incorporated into FWI. Such coherency was earlier used in phase-weighted stacking (PWS) to enhance weak but coherent signals. The benefit of this amplitude-unbiased misfit function is that it can extract information uniformly for all seismic signals (surface waves, reflections, and scattered waves). Using the adjoint-state method, we show how to calculate the gradient of this new misfit function. We validate the robustness of the new approach using checkerboard tests and synthetic data contaminated by random noise. We then apply the new FWI approach to a field dataset acquired at an archaeological site located in Ostia, Italy. The goal of this survey was to map the unexcavated archaeological remains with high-resolution. We identify a known tumulus in the FWI results. The instantaneous-phase coherency FWI results also establish that the shallow subsurface under the survey lines is quite heterogeneous. The instantaneous-phase coherency FWI of near-surface data can be a promising tool to image shallow small-scale objects buried under shallow soil covers, as found at archaeological sites.

Seismic incoherent noise and waves scattered from objects in the crossline directions can cause 2D elastic full-waveform inversion (FWI) to produce artifacts in the resulting 2D models. We develop a complete workflow that can determine subsurface S-wave velocity (VS) models inverted from 2D near-surface seismic data more stably. We make use of a combination of supervirtual interferometry and a matched filter to accurately retrieve dominant surface waves from the field data, whereas the incoherent noise and 3D scattering events are significantly suppressed. The subsurface structures obtained from inverting the retrieved data can be interpreted together with the sections resulting from FWI of the original data to mitigate the potential misinterpretation of artifacts. Our results demonstrate that it is possible to invert 2D near-surface seismic data even when the data quality is lowered by the presence of strong noise and 3D scattered events caused by objects located in the crossline direction.

Geomechanical modelling is generally used to simulate the nucleation of induce d earthquakes in, for instance the Groningen gas field. We apply quasi static simulation to investigate the stress changes from gas production. When a fault reaches a critical state, dynamic simulation provides information on the dynamic rupture during ea rthqu ake nucleation and the resulting wavefield . With the use of geomechanical modelling, it is possible to investigate the effects of the model parameters, e.g., depletion pattern and friction parameters. I n the modelling, the dynamic rupture at a finite fault is simulated both in space and time. The generated seismic wavefield from such a finite source is considered to be more realistic than the resulting wavefield from a point source. T he latter is often assumed in previous studies on the inversion of in duced earthquake data in the Groningen area. To link the wavefield generated by a geomechanically simulated finite source to the field seismic data for an earlier earthquake, we apply the same full moment tensor inversion to the waveform of a finite and of a point source . The inverted moment tensor from the field seismic observation provides a constraint to our geomechanical simulation. This allows us to perform a more realistic simulation of an induced earthquake.