Marmousi Model
Introduction
The AGL Elastic Marmousi model is the most recent generation of the classic Marmousi seismic benchmark. In this tutorial we will cover how one can use this model within Salvus. In particular, we will cover:
- How to read a SEG-Y file using Obspy
- How to read in a rectilinear model using
xarray - How to interpolate a 2-D regularly-gridded model onto the mesh
- How to generate a synthetic seismic shotgather using the model
Environment setup
As always, our first job is to import the relevant Python packages into our notebook.
%matplotlib inline
%config Completer.use_jedi = False
import numpy as np
import os
import pathlib
import obspy
import requests
import time
import xarray as xr
import salvus.namespace as sn
from salvus.namespace import simple_config as sc
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
SIMULATION_MAX_FREQUENCY = 5.0Model building
Information about the Marmousi model we’ll be using can be found here. The model is distributed by the Allied Geophysical Laboratory at the University of Houston as three gzipped SEG-Y files, one each for P-wave velocity, S-wave velocity, and density. The next cell downloads them for you (~5 MB in total).
Marmousi-II was created by Gary S. Martin, Robert Wiley, and Kurt J. Marfurt. Copyright 2004 Allied Geophysical Laboratory, University of Houston, TX, USA. The SEG model page documents its release under CC BY 4.0. The model is provided without warranties; see the license terms.
This tutorial downloads the model from AGL, subsamples it by a factor of ten in each spatial direction, and converts its values to SI units. When sharing the model or model plots, retain the attribution and license information and identify any changes.
Before downloading, check the AGL download page, data access terms, and accompanying notices. If the license for these files is unclear, contact AGL before using or redistributing them.
# Download the Marmousi-II model (~5 MB).
model_directory = pathlib.Path("marmousi-ii")
model_directory.mkdir(exist_ok=True)for name in ["vp", "vs", "density"]:
filename = model_directory / f"{name}_marmousi-ii.segy.gz"
if not filename.exists():
url = f"https://www.agl.uh.edu/downloads/{filename.name}"
response = requests.get(url)
response.raise_for_status()
filename.write_bytes(response.content)Reading the model
We can now start building the model and prepare meshing it in Salvus.
Most material models in seismic exploration are discretized on a regular grid. As models can be stored in a wide variety of formats such as SEG-Y, or RSF, we use xarray as an intermediary to encapsulate both the model parameters and geometry. More info on xarray, including extensive documentation, can be found here.
In general, if you are working with a regularly gridded model in either two or three dimensions, all you need to do to get your model into a form Salvus Mesh can understand is an xarray dataset with the following information
- The coordinate values in meters (named ‘x’, ‘y’, and optionally ‘z’)
- The material parameters in SI units, with names following this parameter naming convention
def read_marmousi(
directory: pathlib.Path, parameters: list[str], decimate: int = 10
) -> xr.Dataset:
"""
Read the Marmousi-II model into an xarray dataset in SI units.
The model is defined on a regular grid with 1.25 m spacing. We decimate
it to speed up reading.
"""
filenames = {"VP": "vp", "VS": "vs", "RHO": "density"}
data = {}
for par in parameters:
stream = obspy.read(
str(directory / f"{filenames[par]}_marmousi-ii.segy.gz"),
format="SEGY",
)
# Each trace is a vertical profile. Velocities are stored in km/s and
# density in g/cm^3, so we convert to m/s and kg/m^3.
data[par] = 1000.0 * np.stack(
[tr.data[::decimate] for tr in stream[::decimate]]
)
# Place the origin in the bottom left corner of the model.
spacing = 1.25 * decimate
nx, ny = data[parameters[0]].shape
model = xr.Dataset(
coords={
"x": np.arange(nx) * spacing,
"y": np.arange(ny - 1, -1, -1) * spacing,
}
)
for par, array in data.items():
model[par] = (("x", "y"), array)
return model# Read and plot marmousi model.
marmousi = read_marmousi(model_directory, ["VP", "RHO", "VS"])A nice feature of xarray.Dataset object is the ability to quickly plot the contained data. Let’s do this below as a way to quickly QC that the model reading went correctly.
marmousi["VP"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))
marmousi["VS"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))
marmousi["RHO"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))<matplotlib.collections.QuadMesh at 0x76aa789d5210>
Note that the top layer models a fluid ocean. In order to accurately incorporate the interface of the ocean bottom into the spectral-element mesh. We will first cut out this part from the model and later on add it back as an ocean layer. This will ensure that the interface conditions between solid and fluid media are correctly imposed and that the model is indeed discontinuous across this interface.
marmousi_elastic = marmousi.where(marmousi["y"] <= 3500 - 37 * 12.5, drop=True)
marmousi_elastic["VP"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))
marmousi_elastic["VS"].transpose("y", "x").plot(aspect="auto", figsize=(15, 5))
marmousi_elastic["RHO"].transpose("y", "x").plot(
aspect="auto", figsize=(15, 5)
)<matplotlib.collections.QuadMesh at 0x76aa76f94d90>
Yup, all looks good!
Initializing a project
We’ll use SalvusProject to manage our workflows. First, we initialize the project (or load it if it is already on disk). Since our simulation domain is defined by the extent of our Marmousi model, we’ll use the from_volume_model constructor to initialize our project.
The ocean layer will be added at the time of meshing, and we will explicitly define the material properties of the ocean later.
project_folder = pathlib.Path("project")
# Initialize a volumetric model from an xarray dataset.
mt = sn.model.volume.cartesian.GenericModel(
name="marmousi_elastic",
data=marmousi_elastic,
)
# Create a project if it does not yet exist or load it from disk.
p = sn.Project.from_volume_model(
path=project_folder, volume_model=mt, load_if_exists=True
)
Up to now, the domain only contains the solid parts, hence its vertical extent is only 3037.5 m.
p.viz.nb.domain()
Next, we also add a purely acoustic version of the Marmousi model to the project. Again, we explicitly model the upper part as a fluid ocean layer to make sure the velocity model is not smeared along the ocean bottom interface.
marmousi = read_marmousi(model_directory, ["VP", "RHO"])
marmousi_acoustic = marmousi.where(
marmousi["y"] <= 3500 - 37 * 12.5, drop=True
)
p.add_to_project(
sn.model.volume.cartesian.GenericModel(
name="marmousi_acoustic",
data=marmousi_acoustic,
)
)Setting up an event
Events in SalvusProject are containers to store waveform data. They are associated with sources and receivers.
Sources
Here we’ll choose a scalar source type (as we’re working in an acoustic medium), and use the side-set interface to allow for easy and intuitive source placement. Note that, since this is a purely acoustic simulation, if we chose an elastic source type (either Vector or MomentTensor) the simulation would give an error explaining that the given source type could not be attached to any element.
Note that at this stage, we only include the spatial characterization of the source. The source time function will enter later, once we specify the frequency of the simulation.
source = sn.simple_config.source.cartesian.ScalarPoint2D(
x=8500.0, y=3490.0, f=1
)Receivers
Setting up receivers is similar to setting up sources. For this example we want to use a regularly-spaced grid of receivers which span the length of the domain. For this, we can use the SideSetHorizontalPointCollection2D receiver type. This allows us to choose a side-set (here, "y1"), and to place an array of receivers a certain distance from this side-set. In this case we choose to place the receivers righ at the ocean bottom. We also need to specify an array of dynamic fields to record. Here we choose pressure, which is equal to "phi_tt".
receivers = (
sn.simple_config.receiver.cartesian.SideSetHorizontalPointCollection2D(
x=np.linspace(0.0, 17000.0, 1000),
offset=0.0,
station_code="xx",
side_set_name="y1",
fields=["phi_tt"],
)
)Because both the source and the receivers are outside of the elastic domain, we skip the consistency check in the following, which would otherwise complain the vertical location.
p.events.add(
sn.EventCollection.from_sources(sources=source, receivers=receivers),
check_consistency_with_domain=False,
)Configuring a simulation
To set up our simulation, we need to specify a few more parameters. The first thing we’ll do is set the duration of the simulation, by specifying end_time_in_seconds.
We could additionally provide parameters like time_step_in_seconds and start_time_in_seconds. If they are not set, Salvus will choose sensible defaults for all three of these parameters. The start and end times will default to the values needed to smoothly simulate the entire source wavelet, and the time step will be computed inside Salvus as per the CFL criterion. Here we set the end time to 8 seconds, and keep the rest of the parameters at their default values.
wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=8.0)We now have enough information to describe an EventConfiguration, which brings together the spatial and temporal aspects of our simulation.
We will be using a Ricker wavelet for this simulation, with a center frequency of half the maximum frequency we’ve set for the mesh.
wavelet = sn.simple_config.stf.Ricker(
center_frequency=SIMULATION_MAX_FREQUENCY / 2.0
)
ec = sn.EventConfiguration(
waveform_simulation_configuration=wsc, wavelet=wavelet
)Now we prepare for adding the ocean layer. The ocean bottom is flat, so the bathymetry model is just a straight line. The width of the ocean column corresponds exactly to the part we cut out at the beginnning of the tutorial.
x = np.linspace(-1e6, 1e6, 10)
y = 37 * 12.5 * np.ones_like(x)
ds = xr.Dataset(
data_vars={"dem": (["x"], y, {"reference_elevation": 0.0})},
coords={"x": x},
)
ocean_layer = sn.bathymetry.cartesian.OceanLayer(
name="ocean",
data=ds,
reference_elevation=0.0,
ocean_surface_datum=3500.0,
ocean_layer_vp=1500.0,
ocean_layer_density=1000.0,
)
p.add_to_project(
ocean_layer,
overwrite=True,
)We’ll also need to specify the absorbing boundary parameters we set above to mock an infinite domain extending from the other side-sets.
# Define coupled Clayton-Enqist / Kosloff
# damping boundaries at the the remaining side-sets.
abc = sn.AbsorbingBoundaryParameters(
reference_velocity=2000.0,
number_of_wavelengths=3.5,
reference_frequency=wavelet.center_frequency,
)With all of our components defined, we can now instantiate a concrete simulation configuration…
Elastic
p.add_to_project(
sn.SimulationConfiguration(
name="marmousi_elastic",
event_configuration=ec,
absorbing_boundaries=abc,
elements_per_wavelength=1.0,
tensor_order=4,
max_frequency_in_hertz=float(SIMULATION_MAX_FREQUENCY),
model_configuration=sn.ModelConfiguration(
background_model=None, volume_models="marmousi_elastic"
),
bathymetry_configuration=sn.BathymetryConfiguration("ocean"),
),
overwrite=False,
)…and then view it below with all the events overlain. Hovering over each source or receiver will reveal more details.
p.viz.nb.simulation_setup("marmousi_elastic", events=p.events.list())[2026-10-12 19:51:05,485] INFO: Creating mesh. Hang on.
<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x76aa76af8a90>
Acoustic
p.add_to_project(
sn.SimulationConfiguration(
name="marmousi_acoustic",
event_configuration=ec,
absorbing_boundaries=abc,
elements_per_wavelength=1.0,
tensor_order=4,
max_frequency_in_hertz=float(SIMULATION_MAX_FREQUENCY),
model_configuration=sn.ModelConfiguration(
background_model=None, volume_models="marmousi_acoustic"
),
bathymetry_configuration=sn.BathymetryConfiguration("ocean"),
),
overwrite=False,
)p.viz.nb.simulation_setup("marmousi_acoustic", events=p.events.list())[2026-10-12 19:51:15,297] INFO: Creating mesh. Hang on.
<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x76aa7f599e50>
Next, we finally launch the simulation in the standard way as below. Remember: if the site you pass is remote, all the shuttling of data to/from the remote machine will be handled for you. Note: since we are using the SalvusProject framework, the simulations will only be submitted if they have not already been. This means that if you run the notebook multiple times, no redundant computations will take place.
The commented lines with the optional argument extra_output_configuration show how to specify volumetric (2- or 3-D) output for wavefield visualization and analysis. We don’t output these by default because the file sizes involved can be quite large (in this case, the parameters below will produce 375 MB of output). If your machine has the capacity, feel free to comment the following in to generate time-varying volumetric output which can be visualized in Paraview.
for sim in ["marmousi_elastic", "marmousi_acoustic"]:
p.simulations.launch(
events=p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=4,
simulation_configuration=sim,
# # uncomment for volume output
# extra_output_configuration={
# "volume_data": {
# "sampling_interval_in_time_steps": 100,
# "fields": ["phi_tt", "stress"],
# },
# }
)[2026-10-12 19:51:18,806] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2610121951039800_6186f3aba3@local
[2026-10-12 19:51:20,819] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2610121951961868_e6a9e7d4ce@local
The previous cell launched our simulations and immediately returned, i.e. it did not wait for the simulations to complete. This is by design, and allows one the flexibility to shutdown and restart the notebook without loosing one’s progress. In this case, however, we do want to wait until the simulations are complete before moving on.
To ensure that we block, we pass an optional arguemnt to query(). This will show update the progress info and wait until we can safely continue with the simulated data. If SalvusProject detects that the simulations have already been run and the results have been stored, the next cell is a no-op.
p.simulations.query(block=True)
True
Visualizing output
Once the simulation has finished, the simulation results will automatically be added to the project.
We can query waveforms for each simulation configuration and event and directly visualize them in the notebook.
p.waveforms.get(data_name="marmousi_elastic", events=p.events.list())[0].plot(
component="A",
receiver_field="phi_tt",
plot_types=["shotgather", "wiggles"],
)
p.waveforms.get(data_name="marmousi_acoustic", events=p.events.list())[0].plot(
component="A",
receiver_field="phi_tt",
plot_types=["shotgather", "wiggles"],
)
Looks pretty cool!