Commercial Application: Seismic Imaging

A seismic survey begins with a physical experiment and ends with a decision made under uncertainty. Between those points, a Linux cluster may read millions or billions of traces, solve the wave equation thousands of times, exchange the edges of three-dimensional domains, preserve selected wavefields, and combine partial results into an image or a model. No one component completes that path. The value appears only when the whole system delivers a scientifically usable result in time for the decision it supports.

This makes seismic imaging a useful closing case for this book. Earlier chapters separated compute nodes, interconnects, storage, facilities, software, scheduling, administration, service operation, and project delivery so that each could be understood. Here we put them back together around one workload. The seismic methods are important, but our subject remains cluster design: how the structure of the problem becomes pressure on memory, communication, storage, and operations.

From a recorded wave to a subsurface model

In an active seismic survey, a controlled source introduces energy into the earth. The energy propagates through rock, changes direction and amplitude at contrasts in physical properties, and is recorded by receivers at known locations. A trace is the sequence of samples recorded by one receiver for one source event. A survey gathers many traces because one source-receiver path cannot describe a three-dimensional subsurface.

The sample values are only part of the record. Processing also needs source and receiver coordinates, sample interval, time or depth reference, component, units, acquisition geometry, and processing history. SEG-Y is a widely used exchange format for storing seismic traces and their headers. The Society of Exploration Geophysicists published SEG-Y Revision 2.1 in 2023 and an extended-textual-header clarification in 2025. A format does not guarantee that metadata are correct, however. Coordinate reference systems, units, scalars, trace identity, and byte order have to be validated before a large run turns a small ingest error into a large wrong result.

Survey geometry multiplies quickly. SEG's open Elastic 2DEW Classic dataset is modest enough to download and inspect: 151 shots, 901 receivers per shot, 2,001 samples per trace at an 8 millisecond interval, occupying 1.1 GB. The much larger SEAM subsalt project reports approximately 65,000 acoustic shot records with as many as 450,000 traces in a shot record. Those figures should not be converted into bytes without the sample format and actual record population, but they show why geometry, samples, and repeated processing matter more than one headline file size.

The recorded traces do not form a picture of the earth. They are observations made at the surface or in wells. Imaging moves recorded energy to plausible subsurface locations using a model of how waves propagate. A velocity model supplies wave speed as a function of position; richer models may include density, anisotropy, attenuation, and elastic parameters. If this model is wrong, a mathematically correct migration can place a reflector in the wrong position.

Migration uses the recorded data and a propagation model to construct an image of subsurface reflectivity. Inversion works in the other direction: it changes model parameters so that simulated data better match the observed data. The two are related but not interchangeable. An image helps an interpreter see structure. An inverted model estimates physical parameters needed to propagate waves and, subject to the assumptions and data, may improve the next image.

One production path

A modern seismic workflow is not one executable consuming one input. Data first enter a controlled project area. Ingest records the source, verifies transferred bytes, parses headers, and produces a geometry and quality report. Preprocessing then corrects or suppresses effects that the selected imaging method is not intended to model. Every correction changes the data presented to the solver, so its parameters and software version belong in the run record.

Model building follows. Well information, interpreted horizons, tomography, previous surveys, and lower-cost inversions may all contribute to a starting model. The team chooses a physical formulation and numerical discretization that are adequate for the question. An acoustic isotropic model is cheaper than an elastic anisotropic model, but the cheaper model can produce a confidently computed answer to the wrong physical problem. The correct choice is the least expensive formulation that preserves the phenomena needed for the decision.

Migration produces an image for quality control and interpretation. Full waveform inversion may then refine the model when the data contain usable information and the starting model is close enough for the optimization to proceed. The refined model can feed another migration. Geoscientists compare gathers, residuals, images, models, and known geology rather than accepting a lower objective function as proof that the earth has been recovered.

Only selected products need to live on the fastest filesystem after the campaign. Raw and irreplaceable inputs, approved models, final images, parameters, source revision, environment description, and quality evidence may require longer retention. Scratch wavefields and trial outputs can often expire. Separating active, checkpoint, project, and archive lifecycles prevents an expensive performance tier from becoming an accidental archive.

Wavefield imaging and its computational consequences

Reverse time migration

Reverse time migration (RTM) is useful where complex velocity structure makes simpler one-way propagation approximations inadequate. For each shot, the method propagates a source wavefield forward through the model. It also injects the recorded receiver data and propagates that receiver wavefield backward in time. An imaging condition combines the two wavefields where and when they coincide. In a common form, the image accumulates their zero-lag cross-correlation at every spatial point and time step. McMechan's 1983 reverse-time migration paper established the basic backward-propagation approach.

The two wavefields explain much of RTM's cluster behavior. A time-stepping solver only needs a few adjacent time levels to continue propagating, but the imaging condition needs the source state corresponding to the receiver wavefield's current reverse-time step. Saving every source time step can exceed accelerator memory and place a large burden on storage. Recomputing the source wavefield reduces stored data but repeats floating-point work. Boundary reconstruction saves selected boundary values and performs another propagation. None of these is free; the best choice depends on memory, device throughput, storage bandwidth, numerical formulation, and job failure policy.

Checkpointing makes the trade explicit. Instead of saving every state, the application retains selected states and recomputes the intervals between them during the reverse pass. The Revolve-based checkpointing work by Kukreja and colleagues describes an optimal schedule for trading extra computation against a smaller memory footprint in inversion problems. On a cluster, the schedule also has a placement question. A checkpoint in device memory is fast but scarce; host memory is larger but crosses the device interface; node-local storage can preserve more state; a parallel filesystem can survive node loss but introduces shared traffic.

An RTM image is normally summed over many shots. Individual shots may be independent through most of their computation, but their images eventually have to be reduced or written for later stacking. A design that treats shot parallelism as perfectly independent can therefore move the bottleneck to the final reduction, file creation, or metadata phase. The complete run includes both the long independent portion and the shorter shared portion.

Full waveform inversion

Full waveform inversion (FWI) uses the wave equation inside an optimization loop. Start with a model, m. For each selected source, solve the forward problem to generate synthetic data. Compare those synthetic traces with the observed traces to form a residual and an objective value. Then solve an adjoint problem driven by the residual. Combining the forward and adjoint wavefields produces a gradient that tells the optimizer how a change in the model would change the objective. Apply a bounded update, and repeat.

This adjoint-state construction avoids computing a separate model derivative for every parameter. It does not make inversion cheap. A single iteration still requires forward and adjoint solves for many sources, plus gradient accumulation, regularization, line search or trial evaluations, and quality control. Plessix's adjoint-state review explains the gradient method, while Virieux and Operto's FWI overview places it in exploration geophysics.

A discrete adjoint or dot-product test checks the implementation before an inversion result is trusted. Apply the forward operator to one test vector and compare the resulting inner product with the inner product obtained by applying the discrete adjoint to a second vector. The two values should agree within a tolerance chosen for the numerical precision and problem size. This test does not validate the earth model, but it can expose an inconsistent forward-adjoint pair that would produce the wrong gradient.

FWI is nonlinear. If the starting model predicts arrivals displaced by too much of a cycle, a local optimizer can reduce the mathematical objective while moving toward an unhelpful model. Missing low frequencies, noise, incomplete acquisition, and physics omitted from the simulator can all make the problem harder. A cluster cannot repair those limitations. More compute permits more sources, frequencies, physical detail, and trials, but each addition has to be justified by information in the data and by validation outside the optimized objective.

The forward wavefield needed for the FWI gradient creates the same store-or-recompute decision seen in RTM. Because an inversion campaign repeats this decision across iterations and trial steps, a small improvement in propagator efficiency or checkpoint policy is multiplied many times. Conversely, a quiet numerical or metadata error can contaminate an entire campaign. Scientific tests belong beside performance tests.

Where the workload becomes parallel

Seismic applications expose several kinds of parallelism, but they do not impose the same costs. A shot ensemble can distribute complete source experiments to independent groups. A time-domain propagation partitions one large spatial grid among ranks and exchanges halo cells at every or selected time steps. A frequency-domain method may distribute independent frequencies, yet the solve within one frequency can require a large complex linear system. Calling all three cases “parallel wave propagation” hides the part that determines the interconnect and memory design.

Shot parallelism is attractive because communication is infrequent during propagation. It is not automatically balanced. Shots can have different receiver counts, recording lengths, preprocessing masks, or convergence behavior. Static assignment may leave a group idle while the longest shot finishes. Dynamic assignment improves utilization, but it has to preserve deterministic bookkeeping and avoid turning a manager into the data bottleneck. The original chapter's PVM example used this manager-worker idea on IBM RS/6000 workstations. PVM is now historical, but the scheduling problem remains visible in MPI work queues, job arrays, and workflow systems.

Spatial domain decomposition is more tightly coupled. Divide a three-dimensional grid into subdomains and each rank can update its interior cells. The stencil near a boundary needs values owned by a neighbor, so ranks exchange halo regions. Computation grows with subdomain volume while communication grows with its surface. Strong scaling eventually gives each rank too little interior work to hide the halo exchange. A thin or badly shaped subdomain reaches that point sooner than a compact one.

This is where the topology and communication path from Chapter 3 become part of the application. Rank placement should keep communicating neighbors close. On accelerator nodes, the active path may include device memory, an on-node scale-up link, PCIe, the network adapter, and the scale-out fabric. Device-aware MPI or another direct path can avoid host staging when the platform supports it, but only a trace of the running application can prove which path was used.

Frequency-domain work has a different shape. Transforming time data into frequency components creates outer parallelism across selected frequencies and shots. In a Helmholtz-based formulation, the work within one frequency includes a complex sparse system whose factorization or iterative solution can consume substantial memory and communication. Frequency-extrapolation migration, including the historical case in the next section, instead advances independent frequency components in depth without solving one global Helmholtz system. A fabric chosen for time-domain nearest-neighbor exchange should not be assumed adequate for either pattern without measuring the actual solver and decomposition.

Gradients and images eventually combine. A global reduction after every small unit of work can erase the benefit of otherwise independent shots. Hierarchical reduction, local accumulation, and less frequent synchronization may help, provided the numerical order and precision remain acceptable. The application should record the reduction method because changing operation order can change floating-point results even when the mathematics is equivalent.

From master-mediated I/O to a parallel data path

The first version of this chapter described a finite-difference, frequency-space post-stack migration. Frequencies could be extrapolated independently, so the compute work already had useful parallel structure. The standard implementation nevertheless placed I/O and distribution on one master. It read frequency data, sent assignments to workers, supplied velocity information, collected migrated values, formed the image, and wrote the result. As the worker count increased, the master and its storage path had to serve more consumers.

It is important to correct one explanation in the old text. MPI-1 did not require applications to perform file I/O on a master rank. It did not define a portable parallel-I/O interface at all. Ranks could call operating-system I/O, but coordinating access and obtaining good performance on the shared filesystems of the period were application problems. MPI-IO arrived with MPI-2, and ROMIO supplied a portable implementation. Master-mediated I/O was a common design choice, not an MPI-1 rule.

The revised implementation allowed ranks to read their assigned frequency data and velocity values in parallel, removing the master's read-and-distribute step from the main path. Rank zero still collected depth contributions for imaging and output in the historical design, so the result was not free of serialization. It was enough to improve the measured runs because it removed a costly serialization point that occurred earlier and more often.

The measurements below were reported for the PARAM 10000, a system described in the original chapter as 40 Sun E450 nodes with four 300 MHz processors per compute node and four NFS servers. They are attributed there to S. Phadke, D. Bharadwaj, and S. Yemeni, “Wave equation based migration and modeling algorithms on parallel computers,” Proceedings of SPG, 1998. We have preserved the recorded inputs and elapsed times; we have not independently reproduced the runs. They are historical application measurements, not estimates of current hardware.

The comparison is valuable because the compute kernel was held conceptually steady while the data path changed. The reduction ranged from about 25 to 35 percent in the four recorded cases. It does not establish that parallel I/O always saves 30 percent. It establishes that an application with abundant compute parallelism can still lose a material fraction of its elapsed time to an avoidable serialization point.

MPI-IO gives modern applications a vocabulary for explicit or individual offsets, blocking or nonblocking access, and collective or noncollective coordination. The MPI 4.1 data-access specification notes that collective operations expose the global access pattern and can therefore permit optimizations unavailable to isolated requests. ROMIO commonly implements collective buffering: selected aggregators combine noncontiguous requests into larger operations better suited to the filesystem. The exact strategy depends on implementation and hints, so collective is an opportunity to optimize, not a guarantee.

The storage system underneath still matters. Chapter 4 described Lustre and the path from clients through metadata and object storage services. A seismic application should align its file layout, request size, stripe choice, concurrency, and retention policy with that path. Thousands of ranks each creating small files may exhaust metadata capacity while using little data bandwidth. One shared file with poorly coordinated small writes may create lock contention. A few enormous transfers can deliver high aggregate bandwidth yet stall downstream readers that need a different organization.

Middleware can separate the application's logical arrays and steps from one physical layout. For example, ADIOS2's BP5 engine supports configurable aggregation and subfiles rather than requiring every producer to write an independent file. It is one design option, not a universal replacement for MPI-IO or seismic formats. Raw exchange data may remain in SEG-Y while solver-native arrays, checkpoints, or streamed analysis use another representation.

Measure before changing the format. Darshan can record application-level POSIX and MPI-IO behavior with low overhead, including access sizes, file counts, time, and rank participation. Filesystem and network counters then show what reached the service. The two views answer different questions: the application trace describes what the code requested; system telemetry describes how the shared infrastructure handled those requests under concurrent load.

Performance follows data movement

Many time-domain seismic propagators apply a finite-difference stencil over a large three-dimensional grid for many time steps. Each update reads neighboring values, performs arithmetic, and writes a new value. Once the working set is larger than cache, memory bandwidth and data reuse can limit the rate before the processor reaches its advertised floating-point peak. Increasing stencil order adds arithmetic and neighbor values; whether this helps or hurts elapsed time depends on the accuracy obtained per grid point and time step, not only on operations per second.

Loop order, spatial blocking, vectorization, and NUMA placement determine how much of the node is used. The original chapter made this point on a 20 MHz IBM RISC System/6000. Restructuring its seismic code for instruction scheduling and data reuse reportedly raised sequential performance from 8 to 26 Mflops. Those numbers have no direct bearing on a present processor, but the reasoning survives. The fastest instruction is unhelpful if its operands repeatedly arrive from the slowest level of memory.

Accelerators make the hierarchy deeper. A GPU can update many grid points concurrently, but its high compute rate must be fed from device memory. The domain and required time states have to fit or be partitioned. Source and receiver operations can be sparse and less regular than the interior stencil. Halo regions move between devices, and saved wavefields move to another memory or storage tier. A kernel benchmark that excludes those operations is useful for tuning the kernel but does not predict an RTM or FWI iteration.

Subdomain shape remains important on GPUs. A compact volume gives more interior work per halo cell, while a thin partition increases communication. Yet the most compact mathematical partition may not map cleanly to the number of devices, memory per device, or on-node topology. The application can use different decompositions across nodes and within a node, but that sophistication requires a placement policy and repeatable topology discovery.

Performance portability adds another tradeoff. Hand tuning can make one kernel excellent on one processor generation while increasing the cost of changing physics or moving to another architecture. Domain-specific systems such as Devito express wave equations symbolically and generate optimized loops. Portable programming layers such as RAJA separate some parallel execution choices from the scientific code. Neither removes architecture-specific measurement; they move part of the optimization and maintenance work into a reusable toolchain.

Published results with their boundaries intact

Seismic performance results are easy to misuse because similar labels can hide different equations, grids, precision, hardware, and measured regions. The examples below are not competitors and should not be used to rank CPUs against GPUs. Each one answers a narrower engineering question.

The Devito paper is especially instructive because it reports both numerical context and hardware behavior. The larger domain contained eight times as many points and ran about 9.5 times longer for the stated operator, within roughly 20 percent of direct volume scaling. This is not a current CPU benchmark. It is evidence about one well-generated stencil on one node while large working sets were constrained by memory behavior.

The Gokhberg and Fichtner study demonstrates a different point. Its 3.5 to 4 times improvement belongs to the authors' heterogeneous spectral-element implementation and comparison baseline. Their reported 60 to 80 percent of peak GPU memory throughput helps explain the result more than the word “GPU” alone. The formulation and the real regional tomography case prevent a direct comparison with the Devito finite-difference timings.

The 2024 LLNL/Chevron eSFWI final report provides a useful example of evidence that is promising but incomplete. It reports greater than 60 percent asymptotic weak-scaling efficiency for the tested adjoint kernels. It also says that the planned reservoir-scale Chevron model was not run and that the integrated framework was tested on a synthetic 2-D Marmousi case. Both statements belong beside the scaling result. Omitting the second would change what the first appears to prove.

Running a seismic campaign

A production run is usually a campaign of dependent jobs rather than one scheduler allocation. Ingest and preprocessing create the inputs. Modeling, migration, or inversion fans work across shots and domains. Reductions and quality checks join the branches. A model update or geoscientist decision may launch the next iteration. The scheduler from Chapter 7 has to support both large tightly coupled jobs and arrays or ensembles without allowing one campaign to make the service unusable for everyone else.

Resource requests should describe the stage. A preprocessing task may be storage-bound and need modest compute. A domain-decomposed propagation may need whole nodes, tightly placed accelerators, and a low-latency fabric. A reduction may need less compute but high read bandwidth. Assigning the maximum resource shape to every stage wastes capacity; splitting the workflow too finely can add queue delay, launch overhead, and fragile handoffs. Measure the end-to-end critical path, not just node utilization.

Failure policy follows the cost of recovery. A short shot can simply rerun. A multi-hour domain-decomposed solve may justify checkpoints. An inversion iteration may be restartable by individual source groups if their contributions are named and committed atomically. A partial gradient must never be mistaken for a complete one after a rank or node failure. Completion manifests, checksums, expected contribution counts, and transactional publication of outputs are more reliable than inferring success from the presence of a large file.

Long campaigns also expose silent faults and configuration drift. Record input identities, geometry validation, model version, source revision, compiler and flags, numerical precision, decomposition, rank and thread placement, container or module environment, libraries, scheduler request, node list, accelerator and driver versions, checkpoint policy, and validation result. Chapter 6's software environments make an executable reproducible only when these run-time choices travel with it.

Reproducibility is not the claim that floating-point bits will remain identical on every architecture. It is the ability to reconstruct the experiment, understand material changes, and recover an accepted result within declared tolerances. The Madagascar reproducible-research system described by Fomel and colleagues joins source, data, processing recipes, and figures. A production organization may use different tools, but the same chain of evidence is valuable.

Operations staff need application-level telemetry as well as component health. CPU or GPU utilization alone cannot show whether an iteration is scientifically progressing. Track completed shots, modeled samples or grid cells, objective evaluations, checkpoint time, bytes moved, failed and retried units, reduction time, and quality gates. Then correlate them with node, fabric, and storage counters. This turns “the cluster is slow” into a bounded question about a stage and a resource path.

The commercial measure is often turnaround rather than peak throughput. Finishing an image two days earlier can be valuable; doubling instantaneous filesystem bandwidth for a stage that occupies one percent of the critical path may not be. Cost includes reserved nodes, licensed software, power, fast storage residency, staff time, failed runs, and the delay before interpreters receive a usable product. The optimal cluster configuration can therefore leave some hardware peak unused if it shortens the complete, validated workflow at lower total cost.

Build an application acceptance case

A procurement benchmark should reproduce the decisions the production workflow will make. It needs a dataset the organization can legally provide to candidates, a reference result, a fixed physical formulation, and a run record. It must separate scientific correctness from performance. A fast answer that fails the adjoint test, moves an imaged reflector, or silently drops traces is a failed run.

The open SEAM Elastic 2DEW Classic data are useful for an accessible control case. Their 151 shots and documented geometry can test ingest, trace bookkeeping, preprocessing, shot dispatch, deterministic reduction, restart, and provenance. At 1.1 GB, however, the dataset fits in memory on an ordinary modern node. It should not be used to claim parallel-filesystem scale. Storage and multi-node domain tests need a synthetic or licensed three-dimensional case whose working set exceeds aggregate cache and whose redistribution pattern resembles production.

Run the correctness control on one node or one accelerator before scaling it. Then hold the problem fixed and increase resources to measure strong scaling. Use weak scaling only with a clearly defined rule for increasing the model or shot workload. Repeat enough times to expose variability, state whether caches are warm, and report failed attempts. A single best run describes an opportunity; an operating commitment needs a distribution.

Component tests remain necessary. Memory bandwidth, GPU communication, MPI latency, and filesystem tests help locate a failure and compare a delivered subsystem with its specification. They do not replace the application case. Chapter 10's requirements-to-acceptance method should connect each seismic deadline and risk to both component evidence and an end-to-end result.

The cluster seen through one application

Seismic imaging makes the architecture in Chapter 1 concrete. Shot ensembles are throughput work until a single source requires several tightly coupled nodes. A large RTM or FWI solve is a capability workload inside a longer throughput campaign. The same project can therefore need broad queue throughput, tightly coupled partitions, and substantial shared storage.

Chapter 2's node balance appears in every stencil update: processor or accelerator throughput, memory bandwidth, memory capacity, NUMA placement, and device topology have to support the chosen domain. Chapter 3's fabric carries halos and reductions. Chapter 4's storage path carries traces, models, checkpoints, and intermediate arrays. Chapter 5's power and cooling limits determine how many high-density nodes can sustain this work rather than merely start it.

The software and operating layers are equally visible. Chapter 6 supplies a controlled numerical environment. Chapter 7 places groups and schedules dependent stages. Chapter 8 keeps nodes, identities, time, telemetry, and configuration consistent during long runs. Chapter 9 turns the machine into a service with data policy, support, allocation, incident response, and lifecycle. Chapter 10 connects the commercial deadline to a tested design.

The historical PARAM 10000 result still belongs in this modern account because it demonstrates a durable method. The application was not made faster by adding processors or replacing the migration equation. Its designers found a serialized path, changed ownership of the data movement, and measured elapsed time again. Modern systems offer much faster devices and richer software, but they reward the same discipline: understand the complete path, preserve the scientific boundary, and optimize the part that delays the accepted result.

References and further reading

1. J. Claerbout, Imaging the Earth's Interior, Blackwell Scientific Publications, 1985.
2. G. A. McMechan, “Migration by extrapolation of time-dependent boundary values,” Geophysical Prospecting, 1983.
3. W. W. Symes and J. J. Carazzone, “Velocity inversion by differential semblance optimization,” Geophysics, 1991.
4. S. Phadke, D. Bharadwaj, and S. Yemeni, “Wave equation based migration and modeling algorithms on parallel computers,” Proceedings of SPG, 1998. Citation and measurements retained from the original chapter.
5. R.-E. Plessix, “A review of the adjoint-state method for computing the gradient of a functional with geophysical applications,” Geophysical Journal International, 2006.
6. J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics, 2009.
7. S. Fomel et al., “Madagascar: open-source software project for multidimensional data analysis and reproducible computational experiments,” Journal of Open Research Software, 2013.
8. A. Gokhberg and A. Fichtner, “Full-waveform inversion on heterogeneous HPC systems,” Computers & Geosciences, 2016.
9. N. Kukreja et al., “High-level Python abstractions for optimal checkpointing in inversion problems,” 2018.
10. M. Louboutin et al., “Devito (v3.1.0): an embedded domain-specific language for finite differences and geophysical exploration,” Geoscientific Model Development, 2019.
11. Lawrence Livermore National Laboratory and Chevron, “Elastic Stochastic Full Waveform Inversion (eSFWI), CRADA No. TC02350 Final Report,” 2024.
12. Society of Exploration Geophysicists, SEG technical standards, including SEG-Y Revision 2.1.
13. SEG Advanced Modeling, open seismic model and data sets.
14. MPI Forum, MPI 4.1 data-access routines.
15. Argonne National Laboratory, ROMIO: A High-Performance, Portable MPI-IO Implementation, users guide.