Finite difference visco elastic modelling

<< Click to Display Table of Contents >>

Navigation:  Modelling >

Finite difference visco elastic modelling

 

Description

This module performs 2D wave-equation forward modelling using a finite-difference (FD) time-stepping algorithm. Starting from a subsurface elastic property model, it propagates a synthetic seismic wave field through the model and records the resulting seismograms at receiver locations. The output can then be used for model validation, synthetic dataset generation, or migration velocity testing.

Two physical approximations are available. The Acoustic scheme models only compressional (P-wave) propagation and is faster. The Viscoelastic scheme additionally accounts for shear-wave propagation and anelastic attenuation (quality factors Qp and Qs), giving more realistic amplitude and waveform behaviour at the cost of higher computation time. Two source-receiver geometries are supported: standard Common Shot acquisition (using real survey geometry) and the Exploding Reflector approximation (for zero-offset modelling from a reflectivity image).

Note: This module is deprecated. It is retained for compatibility with existing processing flows. New projects should use the current finite-difference modelling workflow available in g-Platform.

 

Input data

GatherItem velocity P-wavs

The compressional (P-wave) velocity model of the subsurface, supplied as a 2D gather in which each trace represents a vertical velocity column. This model is mandatory for both Acoustic and Viscoelastic schemes and defines the lateral and vertical extent of the simulation grid. Velocity values should be in metres per second. The grid spacing (distance between traces and sample interval) must be uniform throughout the model gather.

GatherItem velocity S-wavs

The shear (S-wave) velocity model, in the same 2D gather format as the P-wave model. This input is required only when the Viscoelastic calculation scheme is selected. When running in Acoustic mode, this input may be left unconnected. S-wave velocities control the propagation speed of shear waves and together with the P-wave velocities determine the elastic impedance contrasts at each reflector.

GatherItem Ro

The density model (Greek letter rho, "Ro") of the subsurface in kg/m³, provided as a 2D gather matching the spatial extent of the P-wave velocity model. Density is required for both Acoustic and Viscoelastic schemes because it enters directly into the wave equation through the impedance relationship. Physically accurate density contrasts at layer boundaries produce realistic reflection amplitudes in the synthetic seismograms.

GatherItem relax p-wavs

The P-wave quality factor (Qp) model, specifying the anelastic attenuation of compressional energy at each point in the model. Provided as a 2D gather in the same grid as the velocity models. Low Q values (e.g., 10&ndash;30) cause strong absorption and rapid amplitude decay with depth; high Q values (e.g., 100 and above) approach lossless acoustic propagation. This input is optional; if not connected, an internal default Q value is used. It is most important for the Viscoelastic scheme when realistic attenuation behaviour is required.

GatherItem relax s-wavs

The S-wave quality factor (Qs) model, specifying anelastic attenuation of shear energy. Provided as a 2D gather matching the spatial grid. As with the Qp model, this input is optional; if not connected, an internal default is applied. Qs is only physically meaningful in Viscoelastic mode. In many crustal rocks Qs is approximately half of Qp, but both should be set from well log or laboratory measurements when available.

Type output

Selects the modelling geometry and output type. Two options are available:

Common Shot (default) &mdash; The module simulates one shot gather per source location defined in the connected sorted gather index (GGatherIndexVectorItem). Each source fires the selected wavelet, and the wave field is recorded at all associated receiver positions. The result is a set of pre-stack common-shot gathers written to the output file, suitable for use in migration or amplitude analysis.

Exploding Reflector &mdash; The reflectivity image is used as the initial wave-field condition. Every reflector in the model "explodes" simultaneously at time zero, propagating upward. This produces a zero-offset post-stack synthetic section directly and is much faster than running full common-shot modelling. Connect the reflectivity gather to the "Exploding reflector Gather" input instead of the gather index. Absorbing boundary cells (tapers) at the left and right edges of the model are excluded from the exploding field automatically.

Exploding reflector Gather

A 2D reflectivity gather used as the initial wave-field condition for Exploding Reflector modelling. Each sample in this gather represents the initial amplitude at that spatial grid point. This input is visible and required only when Type output is set to Exploding reflector. The gather must have the same spatial grid dimensions as the P-wave velocity model. Leave unconnected when using Common Shot mode.

Indexes of sorting GGather

The sorted gather index (GGatherIndexVectorItem) that defines the Common Shot acquisition geometry: source positions and their associated receiver arrays. The module iterates over each source in this index and propagates the wave field for that shot. This input is visible and required only when Type output is set to Common Shot. Receiver depths are automatically adjusted to match the topography encoded in the P-wave velocity model.

Parameters

Max time of colculation

The total duration of the simulated wave field, in seconds. Default is 3 s. The module time-steps the FD solver from t = 0 up to this value, recording receiver amplitudes at each output sample. Increase this value to capture deep reflections that arrive late; reduce it to shorten computation time when only shallow targets are of interest. Together with the Sample interval, this parameter determines the number of output samples per output trace. The wavelet length is also automatically set equal to this value in the internal solver.

Sample interval

The time sampling interval of the output seismogram traces, in seconds. Default is 0.002 s (2 ms). This is the sample interval that will appear in the output SEG-Y traces. Note that the solver runs internally at a much finer time step (the "delta t in calculation model" parameter below). The module sub-samples the internal wave field every N internal steps, where N = Sample interval / ModelDeltaT. The Sample interval must be an exact integer multiple of ModelDeltaT; if not, the solver automatically adjusts ModelDeltaT to the nearest valid value.

OutputFileName

The base path and filename for the output synthetic seismogram file(s). The module appends the shot index and a file extension automatically when writing multiple common-shot gathers. Provide a full path to a writable location. The output is written in g-Platform internal GSD format and can subsequently be loaded into the seismic viewer or passed to migration modules.

delta t in calculation model

The internal time step of the finite-difference solver, in seconds. Default is 0.0004 s (0.4 ms). This value must satisfy the Courant&ndash;Friedrichs&ndash;Lewy (CFL) stability condition: it must be small enough relative to the spatial grid spacing and the maximum velocity in the model. A rule of thumb is dt < dx / (V_max * sqrt(2)). Using too large a value will cause the FD solution to become numerically unstable and produce noise or diverge. The Sample interval must be an integer multiple of this value. If the condition is not met exactly, the module automatically rounds ModelDeltaT to the nearest valid value and logs a warning.

WaveParams

A group of parameters describing the source wavelet injected at each shot location. Three sub-parameters are available:

ImpulseType &mdash; selects the analytical wavelet shape. Options include Ricker1, Ricker2, AKB, Berlage, Gaussian, GaussianDeriv, MinPhase, Klauder, Ormsby, Spike, Zero, and Unit. Ricker1 (the second derivative of a Gaussian) is the most commonly used for forward modelling because it produces a compact, zero-phase wavelet with a well-defined peak frequency. Spike is useful for impulse-response testing.

Frequency &mdash; the dominant (peak) frequency of the source wavelet in Hz. Choose a frequency appropriate for the target depth and spatial grid spacing. The spatial grid step must resolve at least 5&ndash;10 grid points per minimum wavelength: dx &le; V_min / (N * f_peak), where N is 5&ndash;10. Higher frequencies improve resolution but require finer spatial and temporal grids and therefore much longer computation times.

WaveLen &mdash; the total time window length of the wavelet, in seconds. In practice the module overrides this with the Max time of calculation value, so the wavelet fills the entire modelled record length. This sub-parameter can generally be left at its default.

Number of grid for absorbing

The width of the absorbing boundary layer (sometimes called a sponge or taper zone) in grid cells, applied to all four sides of the simulation domain. Default is 180 cells. This boundary zone gradually damps outgoing waves so they do not reflect back from the edges of the model and contaminate the synthetic record. A wider taper suppresses reflections more effectively but increases the total model size and therefore memory use and computation time. A minimum of 30&ndash;50 cells is typically needed; 180 cells provides strong suppression for most applications. The velocity and density model gathers are automatically padded by this number of cells on all sides before the simulation begins.

Damping constant of taper

Controls the strength of the amplitude damping applied within the absorbing boundary layer. Default is 0.15. At each time step, particle velocities inside the taper zone are multiplied by a cosine-squared taper that ramps from this factor at the outer edge to 1.0 at the interior boundary. Larger values damp more aggressively and absorb energy faster, but if the value is too large relative to the taper width, impedance contrast between the model interior and the outer taper can generate artificial reflections. The default value of 0.15 works well for typical taper widths of 100&ndash;200 cells.

Type upper board

Defines the boundary condition applied at the top of the simulation domain. Two options are available:

Free (default) &mdash; a free-surface boundary condition is applied at the top of the model. This causes upward-propagating waves to reflect back down from the surface, generating surface multiples and ground roll in the synthetic data. Use this option when you want realistic surface reflections and multiples in the output, for example when modelling land or shallow marine acquisition.

Taper &mdash; an absorbing taper is also applied at the top boundary, suppressing surface reflections and multiples. Use this option when primary reflections are desired without interference from surface-related multiples, for example when generating clean primary synthetic data for AVO analysis or to test migration algorithms without multiple contamination. When Taper is selected, the model is extended upward by NumTapers cells to accommodate the absorbing zone.

Calculation Scheme

Selects the physical wave equation used by the finite-difference kernel. Two options are available:

Acoustic (default) &mdash; runs a 4th-order acoustic finite-difference kernel. Only the P-wave velocity and density models are used. This scheme is computationally efficient and sufficient when shear-wave effects and realistic attenuation are not needed. Suitable for most zero-offset and NMO-corrected synthetic generation tasks.

Viscoelastic &mdash; runs a 4th-order visco-elastic finite-difference kernel. All five model gathers (Vp, Vs, density, Qp, Qs) are required. This scheme propagates both P-waves and S-waves simultaneously and applies frequency-dependent attenuation at each grid point. The output seismograms will contain mode-converted S-wave arrivals, realistic amplitude decay with depth, and correct phase distortions caused by absorption. Use this scheme for full-waveform inversion preparation, elastic AVO modelling, or when the medium has significant shear-wave content.

Output data

GGatherItemOut dVx

The horizontal (x-direction) component of particle velocity recorded at receiver positions, output as a seismic gather. This output captures the horizontal motion component of the simulated wave field. It is primarily useful for converted-wave (PS) analysis and for full-waveform quality control when running the Viscoelastic scheme.

GGatherItemOut dVz

The vertical (z-direction) component of particle velocity recorded at receiver positions, output as a seismic gather. For P-wave dominated propagation in an acoustic or weakly anisotropic medium, this is the primary synthetic seismogram output, equivalent to what a conventional vertical-component geophone would record. Use this gather for comparison with field data, migration testing, or amplitude analysis.

GGatherItemOutPS

A combined or auxiliary output gather containing converted PS-wave energy. This output is most informative when using the Viscoelastic scheme, where mode-converted energy is present in the simulated wave field. It allows separate inspection of shear-wave arrivals without interference from the dominant P-wave events in the dVz output.

Visualisation

Five wave-field snapshot views are updated at regular intervals during the simulation (approximately every 300 time steps) and can be monitored in real time while the module is running:

Snap of P-potential &mdash; the divergence of the particle velocity field, which isolates the P-wave (compressional) potential. Bright areas show where compressional energy is concentrated in the subsurface at the current snapshot time.

Snap of S-potential &mdash; the curl of the particle velocity field, which isolates the S-wave (shear) potential. Most useful when running the Viscoelastic scheme to verify that mode-converted energy is behaving as expected.

Snap of Vx and Snap of Vz &mdash; raw snapshots of the horizontal and vertical particle velocity components throughout the entire model grid, including the taper zones. These are useful for diagnosing numerical stability issues or confirming that the absorbing boundaries are functioning correctly.

VelocitiesWithTapers &mdash; a map of the P-wave velocity model after the absorbing taper cells have been added. Use this view to verify that the model padding is correct and that the taper zone extends as expected around all sides of the original model.