Synthetic Seismic Modelling from Classified Meshes

Scope. This page is the methodological companion to the Seismic Modelling user guide. It documents the physics, the numerical scheme, and the software architecture in enough detail to reproduce or extend the method. The user guide covers the workflow and the dialog; here we explain why the result looks the way it does.

Abstract

We describe a forward operator that turns a classified triangulated surface mesh — an outcrop or interpreted horizon model carrying per-vertex facies and acoustic impedance — into a band-limited synthetic seismic section. Contacts between facies are treated as reflecting elements; their reflectivity follows from the impedance contrast across the contact, and their orientation follows from the in-plane impedance gradient under an explicit two-dimensional extrusion assumption. A Kirchhoff diffraction stack sums the contribution of every reflecting element into each trace using straight-ray traveltimes in a constant background velocity, geometric spreading, and an obliquity weight; each trace is then convolved with a chosen wavelet. The same traversal records, co-located with the modelled reflections, the facies and impedance that produced them, which are exposed as filterable attribute layers on the resulting section. We close with the data model, the colour-mapped panel renderer, and the assumptions and limitations of the approach.

1. Introduction

Tying outcrop interpretation to the seismic scale is a recurring problem in reservoir geology: a contact that is obvious on a cliff face may image strongly, weakly, or only as a diffraction once it is band-limited and viewed under a particular acquisition geometry. Rather than exporting to an external modelling package, VRGS computes a synthetic section in situ from the same classified mesh the interpreter is already working with. The method is deliberately lightweight — a single-scattering (Born/Kirchhoff) approximation with straight rays — so that it runs interactively on a background worker and answers the qualitative question "how would this geometry image?" rather than producing a calibrated amplitude prediction.

2. Method

2.1 Inputs

The operator consumes a triangulated mesh with two per-vertex scalar layers:

Symbol	Quantity	Role
$c_i$	facies / class id at vertex $i$	identifies contacts (where $c$ differs across a facet)
$Z_i$	acoustic impedance at vertex $i$	drives reflectivity and reflector orientation

Impedance is the product of bulk density and P-wave velocity, $Z=\rho V_p$. The layers are produced upstream (classification, painting, or the Add Test Facies + Impedance helper) and are read by name; the operator does not classify.

2.2 Reflector extraction

The mesh is traversed facet by facet. For a facet with vertices $\{i,j,k\}$ we group the vertices by facies id. A facet whose three vertices share a single facies is interior to a unit and contributes nothing. Otherwise the facet straddles a contact and becomes a reflecting element with:

• a centroid $\mathbf{x}=\tfrac{1}{3}(\mathbf{p}_i+\mathbf{p}_j+\mathbf{p}_k)$ in true world coordinates;
• a true 3-D area $a=\tfrac12\lVert(\mathbf{p}_j-\mathbf{p}_i)\times(\mathbf{p}_k-\mathbf{p}_i)\rVert$;
• a reflection coefficient and a normal (below).

2.3 Reflectivity

Each facies group has a mean impedance. For a facet touching two facies the contact is unambiguous; for a facet touching three, the pair with the largest impedance contrast is selected (a deterministic, vertex-order tie-break). With $Z_{\mathrm{hi}}$ and $Z_{\mathrm{lo}}$ the mean impedances on the two sides, the normal-incidence reflection coefficient is

$$R \;=\; \frac{Z_{\mathrm{hi}} - Z_{\mathrm{lo}}}{Z_{\mathrm{hi}} + Z_{\mathrm{lo}}}.$$

This is the standard acoustic, normal-incidence form; angle-dependent (AVO) effects are not modelled.

2.4 The two-dimensional extrusion assumption

The reflecting orientation is not taken from the cliff-face (mesh) normal. Instead it is the impedance gradient across the facet, projected into the acquisition plane and normalised:

$$\hat{\mathbf{n}} \;=\; \frac{P_\Pi(\nabla Z)}{\lVert P_\Pi(\nabla Z)\rVert}, \qquad P_\Pi(\mathbf{v}) = \mathbf{v} - (\mathbf{v}\cdot\hat{\mathbf{m}})\,\hat{\mathbf{m}},$$

where $\nabla Z$ is the constant gradient of the linearly-interpolated impedance over the triangle, $\hat{\mathbf{m}}$ is the horizontal normal of the vertical acquisition plane, and $P_\Pi$ removes the out-of-plane component. The gradient is undefined (returns $\mathbf 0$) for a degenerate triangle or vanishing contrast.

Geologically this asserts that each contact extends perpendicular to the section — a 2-D structure swept along strike. It is what makes a horizontal bed reflect strongly to a source directly above it (its normal points up, into the plane), while a vertical fault is weak at zero offset (its normal is near- horizontal) yet still diffracts off its truncated tips. Where measured per-vertex dip vectors are available they can be substituted for the computed gradient, turning the assumption into data.

2.5 The Kirchhoff diffraction stack

Every reflecting element is treated as a Huygens secondary source. For a source–receiver pair $(\mathbf s,\mathbf g)$ defining trace $T$, and an element at $\mathbf x$ with reflectivity $R$, area $a$, and unit normal $\hat{\mathbf n}$, let

$$r_s=\lVert\mathbf x-\mathbf s\rVert,\qquad r_g=\lVert\mathbf x-\mathbf g\rVert.$$

The element deposits an impulse of amplitude

$$A(\mathbf x) \;=\; \frac{R\,a\,\mathcal{O}}{r_s\,r_g}, \qquad \mathcal{O} \;=\; \tfrac12\!\left(\bigl|\hat{\mathbf d}_s\!\cdot\hat{\mathbf n}\bigr| + \bigl|\hat{\mathbf d}_g\!\cdot\hat{\mathbf n}\bigr|\right),$$

where $\hat{\mathbf d}_s,\hat{\mathbf d}_g$ are the unit ray directions from the element to source and receiver. The factor $1/(r_s r_g)$ is the product of two 3-D spherical-spreading terms; the obliquity $\mathcal O$ down-weights elements seen at grazing incidence and reduces to unity when the normal is undefined.

Traveltimes use straight rays in a constant background velocity $V_0$:

$$t(\mathbf x) = \frac{r_s + r_g}{V_0}\quad(\text{time axis}), \qquad d(\mathbf x) = \frac{r_s + r_g}{2}\quad(\text{depth axis}).$$

The impulse is laid into the trace by linear interpolation between the two nearest samples at $\tau=t/\Delta t$ (or $d/\Delta t$). Summing over all elements gives the reflectivity series

$$u_T(\tau) \;=\; \sum_{\mathbf x} A(\mathbf x)\,\Lambda\!\left(\tau - \tfrac{t(\mathbf x)}{\Delta t}\right),$$

with $\Lambda$ the unit triangular (linear-interpolation) kernel. This is a single-scattering approximation: multiples, transmission losses, mode conversion, and intrinsic attenuation are not represented.

2.6 Wavelet and band-limitation

Each finished trace is convolved once with a source wavelet $w$ to produce the displayed, band-limited trace $s_T = u_T * w$. Two wavelets are available:

• a Ricker wavelet of peak frequency $f$, $w(t) = \bigl(1 - 2\pi^2 f^2 t^2\bigr)\,e^{-\pi^2 f^2 t^2},$ windowed to $\pm$ one period and peak-normalised;
• a Klauder wavelet (the autocorrelation of a linear sweep $f_1\!\to\!f_2$ of length $L$), appropriate for modelling Vibroseis data.

Both kernels are symmetric and odd-length, centred on their midpoint. The wavelet bandwidth controls reflection thickness; the vertical resolution is governed by the classic tuning limit $\lambda/4 = V_0/(4f)$.

2.7 Acquisition geometry

The default geometry is a zero-offset line of co-located stations hung above the model, auto-centred on the mesh footprint and scaled to it: the line runs along the longer horizontal axis, the record length $N\Delta t$ is sized so the panel spans from the line down past the base of the model, and the wavelet frequency is set so several wavelengths fit the section (a fixed survey-scale default would place every reflection in a thin band and leave the panel flat). The geometry is fully general, however — fixed offsets or arbitrary source/receiver trajectories are supported, and a vertical fault that is invisible at zero offset can be imaged by shooting offsets that capture its tip diffractions.

The line can also be aligned to a measured structure: choosing a strike shoots the line perpendicular to it (along the dip direction), and the regional dip is then flattened out (§2.8) so dipping beds image as horizontal bands. Both can be entered directly or taken from a nearby orientation measurement (whose dip-azimuth fixes the strike as $\text{azimuth}-90^\circ$, and hence the line as the dip azimuth).

2.8 Co-located facies and impedance attributes (flattened overlay)

Alongside the amplitude, the forward model records the facies id and impedance of the high-impedance side of each contact and exposes them as two auxiliary panels, co-registered with the reflections, that answer "which contact, of what impedance, produced this event?" directly on the synthetic.

Crucially, these attributes are not placed at the Kirchhoff traveltime used for the amplitude. Doing so would smear each contact's facies along its diffraction hyperbola; a winner-take-all colour fill then tracks the apex locus of that smile, so a flat contact would image as an upward arch instead of a horizontal band. Instead the attributes are projected geometrically onto the section, at the contact's true position.

The section frame. The acquisition line $\hat{\mathbf t}$ is shot ~along the dip direction; the depth axis $\hat{\mathbf q}$ is vertical $(0,0,-1)$; the off-plane direction is $\hat{\mathbf m}=\hat{\mathbf t}\times\hat{\mathbf q}$ (the strike). For a contact element at $\mathbf x$ and a trace with source $\mathbf s$, write $\mathbf d=\mathbf x-\mathbf s$ and resolve

$$z = \mathbf d\cdot\hat{\mathbf q}\ \text{(vertical depth)},\qquad u = \mathbf d\cdot\hat{\mathbf t}\ \text{(along-line)},\qquad p = \mathbf d\cdot\hat{\mathbf m}\ \text{(off-plane)}.$$

Flattening. Each element's depth is reduced by its own along-line dip ramp,

$$z_{\text{flat}} = z - u_{\mathbf x}\,\tan\delta,\qquad u_{\mathbf x} = (\mathbf x - \mathbf c_{\text{line}})\cdot\hat{\mathbf t},$$

where $\delta$ is the regional dip and $\mathbf c_{\text{line}}$ the line centre. Because a bed dipping $\delta$ along the line deepens by $u\tan\delta$, every point of that bed maps to the same flattened depth — a thin horizontal band. The element is written at sample $\tau = 2z_{\text{flat}}/(V_0\,\Delta t)$ (time) or $z_{\text{flat}}/\Delta t$ (depth), onto the trace nearest along the line (every trace within a short window, $\sim$two stations, so a dense contact fills without gaps). When several elements compete for one cell, the one with the smallest off-plane distance $|p|$ wins. The deposits are then dilated to the wavelet's main-lobe half-width; empty cells are left undefined. Doing the flatten per element (not as a per-trace post-shift) keeps the band thin — the along-line window does not smear it.

For $\delta = 0$ this reduces to placing each contact at its true depth on the trace above it — horizontal bands, no arch. With a regional dip, dipping bed-parallel contacts collapse to a constant flattened depth and image as horizontal bands. The amplitude is flattened to match by the corresponding per-trace shift (each trace moved up by $u_{\text{tr}}\tan\delta$), so the overlay stays co-registered with the wiggle (approximately for off-plane or strongly-dipping geometry).

3. Implementation

3.1 Data model

The engine (SeismicModeller) is free of UI, rendering, and mesh types: the caller adapts the mesh into a plain MeshView (vertex positions plus the two attribute spans) and receives a TraceGather of amplitude, facies, and impedance grids. The displayable result, CSeismicSection, stores these as flat arrays laid out [sample · n_traces + trace] so they upload directly as textures, and carries the acquisition-line geometry as a relative-to-eye anchor and three corners for floating-point-safe rendering at large coordinates.

3.2 Rendering

The section is drawn as a single textured vertical panel hung from the acquisition line, with sample 0 at the top and time/depth increasing downward. With a regional dip set the gather has already been flattened (§2.8), so the panel itself stays vertical — the flattening lives in the data, not the panel geometry. The active value grid is uploaded as a single-channel floating-point texture; a fragment shader samples it, applies clip/gain/polarity, maps the result through a 256-entry colour lookup, and writes logarithmic depth so the panel composes correctly with the rest of the scene. For a time axis the vertical extent maps two-way time to depth using $V_0/2$.

3.3 Attribute layers, histograms, and filters

Amplitude, facies, and impedance are exposed as sibling attribute layers reusing the application's standard attribute system. Selecting a layer recolours the panel through that layer's own colour map; each layer carries a histogram and a value filter. Filtering masks the panel — samples outside the filter range (and cells with no recorded attribute) are made transparent — so the interpreter can, for example, isolate the reflections produced by a single facies or by impedance contrasts above a threshold. Undefined cells are encoded as not-a-number and discarded in the shader, which is also how the filter mask is realised.

3.4 Concurrency

The forward model runs on a background worker; the result object is published to the scene immediately and its render path is gated until the worker has finished populating the buffers, so the UI stays responsive while a dense mesh is modelled.

4. Parameter reference (theory view)

Parameter	Symbol	Effect
Background velocity	$V_0$	depth↔time mapping and geometric spreading
Sample interval	$\Delta t$	vertical sampling; with $N$ sets the record length
Number of samples	$N$	trace length / imaged depth
Wavelet frequency	$f$	bandwidth → reflection thickness, $\lambda/4$ resolution
Azimuth, spacing, count, height	—	acquisition line position and sampling
Offset mode / distance	—	zero-offset vs. offset acquisition (fault-tip imaging)

5. Assumptions and limitations

1. Single scattering. No multiples, transmission loss, mode conversion, or attenuation; amplitudes are qualitative, not calibrated.
2. Straight rays, constant $V_0$. No ray bending; traveltimes are exact only for a homogeneous background.
3. Normal-incidence reflectivity. No AVO; $R$ ignores incidence angle.
4. 2-D extrusion. Contacts are assumed to extend perpendicular to the section; genuinely 3-D structure off-plane is not represented unless measured dip vectors are supplied.
5. Surface, not volume. Reflectivity comes from facies contacts on the mesh; the medium between contacts is otherwise non-reflective.
6. Sparse attributes. Facies/impedance exist only where reflections were deposited (co-located with events), by design.

6. Future work

Angle-dependent reflectivity (AVO), curved-ray or eikonal traveltimes for a heterogeneous $V_0(\mathbf x)$, multiple-arrival handling, and substitution of measured dip vectors for the computed impedance gradient are natural extensions, as is migrating the inline trace storage to a sidecar for very large gathers.

References

1. Ricker, N. (1953). The form and laws of propagation of seismic wavelets. Geophysics, 18(1), 10–40.
2. Klauder, J. R., Price, A. C., Darlington, S., & Albersheim, W. J. (1960). The theory and design of chirp radars. Bell System Technical Journal, 39(4).
3. Schneider, W. A. (1978). Integral formulation for migration in two and three dimensions. Geophysics, 43(1), 49–76.
4. Bleistein, N. (1987). On the imaging of reflectors in the earth. Geophysics, 52(7), 931–942.
5. Carcione, J. M., Herman, G. C., & ten Kroode, A. P. E. (2002). Seismic modeling. Geophysics, 67(4), 1304–1325.
6. Sheriff, R. E., & Geldart, L. P. (1995). Exploration Seismology (2nd ed.). Cambridge University Press.

Related

• Seismic Modelling — the user workflow and dialog reference.