ROB VESTRUM, PHD, PGEO

THRUST BELT IMAGING, ROB@TBI.CA

TIM MACARTHUR, MSC, PGEO

THRUST BELT IMAGING, TIM@TBI.CA

https://doi.org/10.66733/ZJYH4579

Abstract

In structurally complex environments like the mountains and foothills of the Andes and Rockies, outcropping rocks and differential weathering generate significant seismic-velocity variation in the near surface. In time processing, we correct for these near-surface fluctuations with weathering statics, minimizing the time delays caused by variability in near-surface velocities due to weathering. These statics calculations are made using the first-break energy on shot records—refraction statics—followed by statics that are calculated with the correlation of seismic reflectors—reflection statics. Correcting for near-surface velocity variation with static shifts assumes vertical raypaths through the weathering layer, and this approximation is violated when we have high-velocity rocks outcropping at the surface.

With depth imaging, we have an opportunity to include the refraction velocities from the first breaks in the PSDM velocity model and raytrace through the weathering layers for a more accurate correction for near-surface velocity variation. The refraction solution is only part of the full weathering correction, which is typically followed by a pass of reflection statics. During our testing, we discovered that the reflection statics from the time processing needed to be replaced with depth-domain reflection statics to develop a complete weathering solution for depth imaging. Before calculating reflection statics, we need to correct for the moveout curvature. An offset-moveout correction that is inherently coupled to the PSDM velocity model gives us a complete depth weathering solution (DWS).

We propose a method where we discard all weathering statics from the time processing and perform all near-surface velocity corrections in the depth domain—including reflection statics. We applied this workflow to 2D and 3D datasets from the foothills of the Andes, which show significant imaging improvements below the mountainous areas of the seismic surveys.

Introduction

Seismic data from onshore geologic settings are shot over rough topography and there are variable velocities just below the acquisition surface. The near surface is the lens through which we image the subsurface. Improved understanding of the lens effects creates an opportunity to more accurately correct for near-surface velocity variation in seismic imaging.

The traditional method for correcting for near-surface weathering is with statics corrections (Lawton 1989, Cox 1999, Yilmaz 2001). Variability in near-surface velocities results in time delays that vary between surface locations of seismic sources and receivers. We represent near-surface velocity variation with static time shifts at each source and receiver location. This practice assumes the raypaths through the near surface are near vertical. This approximation is accurate for a thin weathering layer with a velocity much lower than the subsurface velocity. In cases where high-velocity strata outcrop at surface or if we want to use more of the near-surface velocity information with a thicker low-velocity layer, then we need an alternative method for integrating the near-surface velocities into our seismic-imaging workflow.

Prestack depth migration creates an opportunity to use more of the near-surface velocity information and improve the seismic imaging of subsurface structures.

Statics corrections in time processing

Refraction Statics

The first-arrival energy on shot records provides velocity information for the near surface. As waves propagate along the surface, the traveltime to each receiver location increases with increasing distance from the source. For a range of receivers near the source, the distance between

Figure 1. Wavefront diagram illustrating refracted waves in the two-layer case when v₂ in the lower layer is greater than v₁ in the upper layer. The colour changes along each wavefront indicate whether the portion of the wavefront is downgoing in the upp upper layer (green), downgoing in the lower layer (blue), or the upgoing refracted portion of the wavefront (red). The upper ray indicates the raypath for the direct arrival and the lower ray indicates the refracted arrival. The place where energy along these two raypaths arrive at the same time is called the critical offset, Xc.

source and receiver divided by the time for energy to arrive at the receiver is the velocity of the strata directly below the surface. At longer offsets, the minimum traveltime observed at the receivers is from energy that followed a raypath through higher velocity strata deeper in the subsurface. Figure 1 shows a diagram of a raypath that has followed along a deep refracting layer. At the critical offset, Xc, the raypath through the refracting layer arrives at the same time as the direct arrival along the surface. For offsets greater than Xc, the slope of the first arrivals on the shot record (Figure 2) is the seismic velocity of the refracting layer.

Figure 2. Example of direct-arriving energy and refracted energy on a shot gather. The source location is at the apex of the shot record and the offset increases outward from the apex. The slope of the first-arriving energy along the nearest offsets is the velocity, v₁, in metres of offset per second of recorded time. For offsets past the critical offset, Xc, the slope of the first arrivals is the refractor velocity, v₂.

There is a variety of algorithms used to obtain velocity information from first-break energy (e.g. Lawton, 1989, Becerra et al, 2009). The inherent redundancy of the seismic experiment—with many sources and receivers covering the same area—creates an opportunity to use statistical inversion algorithms to find optimized velocity structures. These inversion methods may assume a certain specified number of layers and invert for layer thickness and velocity. In areas of complex structure, the more robust algorithms are the grid-based first-arrival tomography methods that calculate a grid of the velocities in the subsurface, without the layering assumption.

An example of a velocity model that results from grid-based tomography method is shown in Figure 3. The workflow from here is to pick a base of weathering and then calculate static shifts based on the differences between the model velocity and the replacement velocity. The replacement velocity is usually chosen as a velocity approximately equal to the highest velocity near the base of the refraction model. The application of statics solution minimizes the effects of variable time delays in the near surface by effectively “replacing” the low velocities with the replacement velocity.

Figure 3. 2D velocity model calculated from first-arrival tomography.

Calculating time shifts vertically across the low-velocity layer uses the approximation that the raypaths through this low velocity layer are vertical. Figure 4 shows a ray diagram for a simple horizontal reflector below a low-velocity layer. When the velocity decreases at the boundary between the subsurface and the weathering layer, the rays refract towards the normal vector of the interface. The statics approximation assumes a small difference between a vertical raypath and the refracted raypath, which is more accurate with a relatively thin weathering layer and a greater difference between the subsurface and weathering velocities.

Figure 4. Ray diagram for a horizontal reflector beneath a weathering layer. The source and receiver locations are marked S and R. S’ and R’ are the source and receiver positions projected vertically onto the base of weathering, illustrating the approximation inherent to statics correction.

Reflection statics

The redundancy of the seismic experiment creates an opportunity for a secondary weathering correction calculation called “reflection statics.” The reflection statics are calculated after refraction statics are applied. The goal is to correct for time delays left behind by the refraction-statics algorithm.

Seismic traces are sorted into the common midpoint (CMP) domain and corrected for normal moveout (NMO) to create flat events on the CMP gathers. Cross-correlations between a collection of CMP gathers calculate shifts required to align the traces. Using these calculations, the reflection-statics algorithm finds shifts that have a consistent trend for each source and receiver location. Constraining the shifts to sources and receivers is called the surface-consistent constraint and the resulting statics are then considered time delays localised to each source and receiver location.

Figure 5 shows an example of NMO-corrected CMP gathers before (Figure 5a) and after (Figure 5b) the application of surface-consistent reflection statics. The reflectors on these gathers look much more coherent with only the alignment correction applied from the reflection statics. Summing static-corrected traces into the final stacked section will show higher coherency on the final imaged reflectors.

Figure 5. Common midpoint (CMP) gathers (a) before and (b) after the application of reflection statics.

The full time-weathering solution (TWS) has a combination of refraction statics followed by reflection statics to minimize the effects of variation in the near-surface velocities.

Weathering velocities in the PSDM velocity model

The velocities we calculate from first-arrival tomography have proven useful for the calculation of weathering statics, and we have an opportunity to use more of that velocity information deeper into the model and velocities higher than the weathering velocities used in statics. The intention was to remove the refraction statics from the input shots and correct for near-surface, low-velocity variation by including the velocities in the PSDM velocity model. Simply substituting the refraction velocities for the refraction statics resulted in mixed, if not degraded imaging results.

One of the testing datasets is a 2D line from the Llanos foothills from MacArthur et al. (2022). This line shows clear degradation of imaging with the simple substitution of refraction velocities for refraction statics. Figure 6a shows a PSDM image from the traditional method where the time-processing statics were applied to the input to the migration. When we included the velocities in the PSDM velocity model and removed the refraction statics from the input, the resulting image (Figure 6b) shows degradation in reflector coherency, particularly in the shallow section.

The reflection statics are calculated on seismic data with refraction statics applied. After replacing the refraction statics with refraction velocities, the inaccuracies of the refraction solution have changed—for better or worse. So, we need a new method to calculate reflection statics that are coupled to the new refraction method.

The total weathering solution for time processing is the refraction statics plus the reflection statics, so the next step in replacing the time weathering solution was to also replace the reflection statics using the model-based moveout technique published by Ellison (2018). This data example underscores the need for a reflection-statics calculation coupled to the depth-migration velocity model.

Figure 6. Prestack depth migration of a 2D seismic line with (a) the refraction and reflection statics from time processing applied on input and (b) the refraction velocities in the PSDM velocity model and the reflection statics from time processing applied to the input. The velocity overlay above each section illustrates the differences in velocities with (b) and without (a) the refraction velocities in the model. The boxes indicate the greatest difference between the two seismic images: the degradation of imaging in the shallow section on (b).

Reflection statics in the depth domain

An essential ingredient to the calculation of reflection statics is the moveout calculation. As shown above, in time processing the moveout calculation is from NMO, which uses the NMO equation and moveout velocities to flatten the CMP gathers. The NMO calculation and its moveout velocities use the hyperbolic assumption and are inherently coupled to the time-processing workflow (Ellison, 2018, MacArthur et al, 2022, Beck et al, 2025). The alternative to NMO in depth is model moveout (MMO), which uses the traveltimes from the PSDM velocity model to calculate the moveout correction.

Figure 7 shows the MMO workflow. Prestack migration does two things: (1) moves the reflection energy into the position of the reflectors and (2) corrects for the moveout effect of increasing source-receiver offset. When we calculate MMO, we separate out these two effects and only apply the moveout correction. One could consider MMO an offset-dependent time-to-depth conversion.

Once we apply MMO, we find our CMP gathers in the depth domain, and our goal is to calculate static shifts that we may apply to the traces input to the depth migration, which are in time. A CMP-constant depth-to-time conversion transforms the gathers back into the time domain. From here, the standard reflection-statics algorithm calculates the required surface-consistent reflection statics. The inputs to the depth migration are the shot gathers with MMO applied and the traveltimes calculated from the velocity model. The result is a statics solution that is inherently coupled to the PSDM velocity model. The depth weathering solution (DWS) has the refraction velocities integrated into the PSDM velocity model and MMO reflection statics on the input gathers.

Figure 8 shows PSDM images with the TWS (Figure 8a) as compared to the DWS (Figure 8b). Note that the workflows that create these two images are identical. The only difference between the two velocity models is in the low-velocity layer in the shallow section. The table below describes the differences between the two PSDM images.

TWS DWS
refraction statics from first-arrival tomography velocities from first-arrival tomography
reflection statics from cross-correlation of NMO gathers statics from cross-correlation of MMO gathers

Figure 7. Flow chart for the model moveout (MMO) workflow. Inputs are the PSDM velocity model and the seismic shot gathers from time processing. The output are depth migrated data with MMO reflection statics applied to the input. These reflection statics are inherently coupled to the PSDM velocity model.

Including MMO reflection statics in the DWS solution resolved the shortcomings of the PSDM imaging in the shallow section (Figure 6) and shows improved imaging in the footwall of the thrust over the TWS (Figure 8).

Figure 8. Prestack depth migration of a 2D seismic line (a) with the time weathering solution applied and (b) with the full depth weathering solution applied. The velocity overlay above each section illustrates the differences in near-surface velocities. The boxes highlight the reflectors in the footwall of the overthrust, which show improved reflector coherency and sharper truncations on (b).

3D data example

Optimized for three-dimensional imaging of subsurface targets, 3D seismic surveys on land have lower data density in the near surface than the single-azimuth acquisition of 2D seismic data (Vestrum and Gittins, 2009). With lower density coverage in the near surface, weathering corrections for 3D surveys often rely more heavily on the reflection statics than 2D seismic data. Figure 9 shows comparisons between TWS and DWS depth images from a 3D survey in the foothills of the Andes. This inline is from the central part of the 3D and has strong reflectivity throughout. The small box in this comparison shows improved structural detail and fault-plane imaging on the DWS image (Figure 9b) as compared to the TWS result (Figure 9a). The larger box shows overall improved reflector coherency of subthrust reflectors.

Figure 10 shows an inline near the edge of the survey that struggles with decreased fold and increased structural complexity subthrust. The DWS (Figure 10a) shows more differences compared to the TWS (Figure 10b), with significantly improved reflector coherency in the foreland reflectors in the tall box to the right of each image. The central box in this comparison (Figure 10) shows improved coherency and resolution of subthrust reflectors, similar to the result on the central crossline.

Discussion and conclusions

In seismic imaging, the near surface is the lens through which we image the subsurface. More accurately correcting for these lens effects results in improved seismic imaging in complex-structure areas like the foothills of the Andes mountains.

During the development of this workflow, we found that simply including the first-arrival tomography velocities in the PSDM velocity model did not yield the results that we had hoped. A workflow that includes reflection statics coupled to the PSDM velocity model offers a complete depth weathering solution, or DWS, that is independent of any weathering corrections or statics solutions from the time processing. The results show improvement in seismic imaging on 2D and 3D seismic data with DWS as compared to the traditional method of using weathering solution from the time processing.

Now that we have a method to correct for the effects of near-surface velocities that is independent of the time-processing solution, we have the opportunity for further workflow optimization. It may be helpful to pick first arrivals to longer offsets and use more of the first-arrival tomography solution in the PSDM velocity model. We may also experiment with more advanced techniques. Diving-wave full-waveform inversion may leverage more of the early-arrival energy on shot records to build a subsurface velocity model that has five to ten times more thickness than our current near-surface models from first-arrival tomography.

Figure 9. Inline A near the centre of a 3D PSDM volume with (a) time weathering solution applied and (b) depth weathering solution applied. The velocity overlay near the top of the inline slice shows the shallow velocities in each velocity model. Boxes indicate key differences in imaging between the TWS and DWS versions.

Figure 10. Inline B near the edge of a 3D PSDM volume with (a) time weathering solution applied and (b) depth weathering solution applied. The velocity overlay near the top of the inline slice shows the shallow velocities in each velocity model. Boxes indicate key differences in imaging between the TWS and DWS versions.

References

Becerra, C., Agudelo, W., Guevara, S., 2009. Uncertainty Analysis in Statics Corrections Obtained by Tomographic Inversion: Application in a Mountainous Zone in Catatumbo (Colombia). The Leading Edge, 28(2): 212–2

Beck, J., Krishnasamy, T., Reta-Tang, C., and Ogunsuyi, F., 2025, Rethinking Statics with FWI, GeoConvention, Calgary, Canada.

Cox, M., 1999, Static Corrections for Seismic Reflection Surveys: Society of Exploration Geophysicists, Tulsa, USA.

Ellison, D.K., 2018, Depth imaging with reflection statics derived from model-based moveout, MSc thesis, University of Calgary.

Ellison, D and Innanen K., 2016, Improved resolution in depth imaging through reflection statics corrections derived from model-based moveout: CREWS Research Report, Vol 28.

Lawton, D.C., 1989, Computation of refraction static corrections using first-break traveltime differences, Geophysics, 54(10), 1289-1296.

MacArthur, T., Vestrum, R.W., and Cameron, G.H., 2022, Correcting for near-Surface velocity variation in seismic depth imaging using model moveout (MMO), AAPG ICE 2022, Cartagena, Colombia.

Murphy, G. H., and Gray, S. H., 1999, Manual seismic reflection tomography: Geophysics, 64, 1546-1552.

Vestrum, R.W., and Gittins, J.M., 2009, Technologies from foothills seismic imaging: replacements or complements?: First Break, 27, (2), 61-66.

Yilmaz, O., 2001, Seismic data analysis: processing, inversion and interpretation of seismic data (2nd ed.): Society of Exploration Geophysicists, Tulsa, USA.

About the Authors

Rob Vestrum’s geophysics education includes a BSc from the University of Alberta and a PhD from the University of Calgary. He worked in the foothills team at Shell and the depth-imaging team at Veritas before founding Thrust Belt Imaging (TBI) in 2004. TBI specializes in seismic imaging for structurally complex settings, with an emphasis on thrust belts and anisotropy.

He has taught courses on seismic imaging and anisotropy for more than 30 years at venues worldwide and has served as a sessional instructor for a graduate course at the University of Alberta. Rob also brings extensive hands-on experience in seismic processing across a wide range of basins.

Recognition for Rob’s contributions to geophysics include the John Kendall award for top doctoral thesis in the Faculty of Science at the University of Calgary, the SEG’s J. Clarence Karcher Award for significant contributions to science by a young geophysicist, and the CSEG Roy O. Lindseth Medal.

Tim MacArthur is a senior geophysicist at TBI with over a decade of global experience in seismic imaging in complex geologic settings. His work involves applying data analytics, visualization, and research-driven innovation to tackle exciting imaging cases. His seismic-imaging experience includes 2D and 3D projects from across the globe, working in both Calgary and Bogotá offices. He has led teams, trained colleagues, driven research and development initiatives, and met the expanding challenges of TBI’s clients.

Tim holds a BSc in Geophysics from the University of Calgary (2013) and an MSc in Information Systems, specializing in Data Analytics, from Athabasca University (2025).

Recent work focuses on developing a novel Reverse Time Migration (RTM) algorithm and Full Waveform Inversion (FWI) designed specifically for rough topography, pushing the boundaries of imaging accuracy in these difficult environments.

Tim is passionate about solving seismic imaging problems and creating lasting impact in energy and beyond.