Mondaic Documentation
www.mondaic.com
Mondaic
Version:
No results found.

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
As always, our first job is to import the relevant Python packages into our notebook.
Copy
%matplotlib inline
%config Completer.use_jedi = False

import numpy as np
import os
import shutil
import pathlib
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.0
Information about the Marmousi model we'll be using can be found here. The model is stored as a SEG-Y file, and the next cell will download and unpack the model for you. Note that the model file is relatively large (~150 MB), so make sure you have enough disk space before continuing.
# Download Marmousi model ~150 MB.
target_path = "elastic-marmousi-model.tar.gz"
if not pathlib.Path(target_path).exists():
    url = (
        "https://s3.amazonaws.com/open.source.geoscience/open_data"
        "/elastic-marmousi/elastic-marmousi-model.tar.gz"
    )

    response = requests.get(url, stream=True)
    if response.status_code == 200:
        with open(target_path, "wb") as f:
            f.write(response.raw.read())
        shutil.unpack_archive(target_path)

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
# Read and plot marmousi model.
model_directory = pathlib.Path("./elastic-marmousi-model/model")
marmousi = sn.toolbox.read_elastic_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 0x74fbf0947810>
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 0x74fbef6e3ad0>
Yup, all looks good!
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 = sn.toolbox.read_elastic_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,
    )
)
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,
)
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...
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())
[2025-11-12 21:29:49,305] INFO: Creating mesh. Hang on.

<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x74fbf550a5d0>
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())
[2025-11-12 21:29:53,298] INFO: Creating mesh. Hang on.

<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x74fbeea58450>
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"],
        #     },
        # }
    )
[2025-11-12 21:29:55,167] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2511122129310863_f91f33c6c2@local

[2025-11-12 21:29:56,599] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2511122129616613_108096344b@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
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!
Previous
Heterogeneous models and slice extraction
Next
3D surface topography
We use cookies to analyze our traffic, but won't track you until explicit consent is given. Using our contact form will create relevant cookies and transmit data to third party services to keep track of your request. Read more about how we use personal information in our Privacy Policy.
PAGE CONTENTS